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lAAMTHACT  (MnAiwtnlWMMW 


The  hybrid  stress  method  has  demonstrated  many  improvements  over  conventional  displacement-based  el¬ 
ements.  A  main  detraction  from  the  method,  however,  has  been  the  higher  computational  cost  in  forming 
element  stiffness  coefficients  due  to  matrix  inversions  and  manipulations  as  required  by  the  technique.  By 
utilizing  special  transformations  of  initially  assumed  stress  fields,  a  spanning  set  of  orthonormalized  stress 
modes  can  be  generated  which  simplify  the  matrix  equations  and  allow  explicit  expressions  for  element 
stiffness  coefficients  to  be  derived.  The  developed  methodology  is  demonstrated  using  several  selected  2-D 
quadrilateral  and  3-D  hexahedral  elements. 
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1  Introduction 


Alternative  approaches  to  displacement-based  element  formulations  began  with  the  pioneering  work  of  Pizm 
[1]  who  developed  the  stress-based  approach  for  computing  element  stiffness  matrices  based  on  the  principle 
of  minimum  complementary  energy.  The  development  of  the  hybrid  stress  technique  followed  using  the 
Bellinger- Reissner  functional  [2]  wherein  both  stresses  and  displacements  are  assumed  as  independent  quan¬ 
tities  and  from  which  the  purely  stress-based  approach  can  be  formally  derived.  Current  variational  bases  for 
finite  element  formulations  range  from  single  field  functionals  to  the  most  general  Hu-Washizu  functional  in 
which  displacements,  strains  and  stresses  are  all  assumed  as  independent  quantities  [3,4].  The  use  of  multiple 
independent  field  variables  in  element  formulations  has  created  a  rich  arena  of  theoretical  approaches  with 
which  to  maximize  finite  element  performance  and  has  yielded  elements  with  improved  convergence  behavior 
and  stress  prediction,  avoidance  of  locking  in  constrained  media  problems,  and  the  inherent  capability  to 
represent  traction-free  edge  conditions  and  singular  stress  fields.  However,  this  has  resulted  in  additional  the¬ 
oretical  requirements  in  the  formulation  of  robust  elements  which  have  been  addressed  by  various  researchers. 
These  include  the  requirement  of  element  invariance  under  coordinate  transformation  [5,6],  suppression  of 
spurious  zero  energy  modes  [7],  minimum  e.\pansions  for  the  independent  field  variables[7,8],  and  optimal 
sampling  points  for  stress  recovery  [5].  A  major  detraction  of  element  formulations  based  on  multi-held  func¬ 
tionals  over  displacement-based  elements  has  been  the  computational  cost  associated  with  matrix  inversions 
and  additional  numerical  manipulations  required  in  generating  element  stiffness  matrices.  As  used  herein, 
computational  cost  or  efficiency  is  meant  to  refer  to  the  minimum  number  of  sequential  operations  formally 
required  in  mathematical  statements  and  not  to  the  efficiency  of  specihc  algorithmic  implementations.  The 
present  work  focuses  on  the  hybrid  stress  method  in  which  the  structure  of  the  element  matrices  is  dehned  by 
the  Hellinger-Reissner  functional.  A  novel  procedure  is  detailed  herein  for  minimizing  the  numerical  cost  by 
making  full  use  of  the  freedom  in  selecting  and  manipulating  assumed  stress  helds  which  enable  explicit  forms 
of  element  stiffness  matrices  to  be  derived.  The  developed  procedure  is  based  on  a  detailed  examination  of 
the  complementary  energy  matrix  inherent  to  the  hybrid  technique  and  leads  to  stress  field  transformations 
which  include  an  orthonormalization  approach  hitherto  not  attempted.  The  application  of  the  technique  is 
demonstrated  with  2-D  quadrilateral  and  3-D  hexahedral  elements  incorporating  incompatible  displacement 
fields.  By  simplifying  the  constituent  matrices  involved  in  forming  element  stiffness  coefficients  such  that 
an  explicit  evaluation  can  be  accomplished,  the  hybrid  stress  method  is  shown  to  offer  a  computational 
advantage  over  similar  displacement-based  element  formulations. 


2  Variational  Basis  for  the  Hybrid  Stress  Model 


The  extended  Hellinger-Reissner  functional  may  be  stated  as 

11/1=  f  [(-l/2)(r^ Scr  +  (T^(Lu)i,)]dv  -  f  vi^tds  (1) 

Jv  Js 

where  <r  is  the  assumed  stress  field,  S  is  the  material  compliance  matrix,  u,  and  ua  are  the  assumed  compatible 
and  incompatible  displacement  fields,  L  is  the  differential  operator  relating  struns  to  displacements,  and  t 
are  applied  surface  tractions  over  a  portion  of  the  element  boundary,  s. 

The  assumed  stresses  may  be  represented  by 

<r=P0  (2) 

where  P  is  a  matrix  of  polynomial  terms  and  /?  is  a  vector  of  undetermined  expansion  coefficients.  Each 
independent  stress  mode  is,  therefore,  represented  by  a  column  in  P.  The  displacement  field  is  assumed  over 
the  element  domain  as 

u  =  +  ux  =  Nq  -t-  MX  (3) 

where  N  and  M  are  compatible  and  incompatible  displacement  shape  functions,  respectively,  q  are  nodal 
displacements,  and  A  are  Lagrange  multipliers  which  enforce  internal  constraints.  In  the  form  of  (1),  the 
incompatible  displacements  act  to  variationally  enforce  the  orthogonality  of  the  stresses  to  the  incompatible 
strain  modes  in  a  weak  sense.  Neglecting  applied  tractions  and  substituting  (2)  and  (3)  into  (1)  yields 


H/i  =  /  [(-1/2)13'^  P'^SPf)  +  P'^(LN)q-^  13'^  P'^{LM)X]dv 

Jv 


(4) 
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or 

Hr  5  {-1/2)0'^ H0  +  0^Gq  -1-  0'^RX 

Jfi 

(5) 

where 

.'.j; 

//"=  f^P'^SPdv 

(6) 

G  =  lP^{LN)dv  =  lP'^Bdv 

(7) 

R  =  f^P^iLM)dv  =  lP'^Bdv 

(8) 

Seeking  a  stationary  value  of  the  functional  by  taking  the  first  variation  with  respect  to  and  A  ^  .elds 


H-^{Gq-  R\)  (9) 

and 

R^I3  =  Q  ■  (10) 

By  eliminating  A  and  substituting  the  resulting  expression  for  /?  into  (5),  the  vziriation  with  respect  to  q 
yields  the  element  stiffness  matrix  as 

K-G'^H-^G  (11) 

where 

G  =  [I -R(R^H-^R)-^R^H-^\G  (12) 


3  Computational  Minimization  in  the  Hybrid  Technique 

Since  the  introduction  of  the  hybrid  stress  method,  minimizing  the  computational  cost  associated  with  equar 
tion  (11)  has  been  an  ongoing  concern.  Counterbalanced  with  the  ostensibly  greater  computational  cost  and 
demonstrated  improvement  in  element  behavior  afforded  by  the  two-field  hybrid  stress  technique  has  been 
the  computational  simplicity  of  displacement-based  elements.  Thus,  the  selection  of  element  formulations 
has  been  a  form  of  Occam’s  Razor  in  what  minimum  degree  of  computational  cost  is  required  to  implement 
a  useful,  convergent  element  to  obtain  accurate  solutions  to  practical  problems.  To  make  hybrid  element 
formulations  competitive,  various  approaches  have  been  applied  to  minimize  the  cost  in  evaluating  equation 
(11).  In  formulations  involving  only  compatible  displacements,  it  has  been  noted  that  is  not  required 
separately  but  the  product  H~^G  can  be  obtained  via  equation  solver  techniques  using  equation  (9)  by 
treating  the  columns  of  G  as  multiple  right  hand  sides  which  leads  to  a  significant  reduction  in  computation 
[5].  Other  simplifications  have  been  achieved  by  the  use  of  functionals  such  as  (I)  in  which  incompatible 
displacement  modes  are  used  to  variationally  enforce  equilibrium  or  orthogonality  construnts  [2,8,9].  This 
permits  the  use  of  unconstrained  and,  therefore,  uncoupled  stress  expansions  which  lead  to  block  diagonal 
representations  of  the  ^-matrix  in  which  the  calculation  of  can  be  limited  to  inversions  of  the  submatrix 
blocks.  Recently,  a  theoretical  basis  has  been  developed  for  making  admissible  variations  of  terms  in  the 
complementary  energy  matrix  which  permit  simplifying  approximations  to  be  made  based  on  stability  and 
convergence  considerations  to  minimize  the  cost  of  computing  [10].  Nevertheless,  the  various  treat¬ 
ments  of  equation  (11)  still  represent  a  significant  computational  cost  in  terms  of  numerical  integrations  and 
manipulations.  The  aim  of  the  present  effort  is  to  strictly  adhere  to  the  variational  constraint  expressed  by 
equation  (11)  by  simplifying  the  fundamental  mathematical  statements  through  formal  procedures  without 
introducing  additional  assumptions  or  approximations. 

4  Assumed  Stress  Expansions 

The  conventional  procedure  in  the  hybrid  stress  method  is  to  define  stress  expansions  in  the  natural  or 
mapped  coordinate  system.  This  definition  has  been  used  to  develop  rational  procedures  for  eliminating 
spurious  kinematic  modes  and  maintaining  element  invariance  while  keeping  the  number  of  independent  stress 
modes  to  a  minimum.  However,  one  drawback  of  this  approach  has  been  the  contravariant  transformation 
required  to  express  stresses  in  physical  coordinates.  This  transformation  will,  in  general,  cause  coupling 
between  the  constant  and  higher-order  stress  terms  and  the  element  will  subsequently  fail  the  patch  test. 
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A  solution  to  this  problem  has  involved  the  use  of  an  approximation  to  the  Jacobian  by  using  its  constant 
value  at  the  element  centroid.  This  approximation  -  or  ’variational  crime’  •  is  reasonably  accurate  for 
constant  strain  elements  with  linear  interpolation  functions  but  should  be  expected  to  demonstrate  increasing 
inaccuracy  with  general  element  distortion  and  element  order  and  be  a  limiting  factor  in  coarse  mesh  accuracy. 
In  the  present  development  stresses  will  be  assumed  in  both  physical  and  njttural  coordinates  to  demonstrate 
a  procedure  for  minimizing  computational  cost  in  all  basic  formulations  using  the  hybrid  technique.  In  using 
stresses  defined  in  physical  coordinates  all  ad  hoc  simplifications  will  be  avoided  and  ’exact’  expressions 
will  be  derived  within  the  formal  approximation  framework  based  on  the  order  of  element  interpolants 
and  variational  constraints.  This  is  done  in  anticipation  of  future  study  in  developing  explicit  higher- 
order  element  matrices  where  the  degradation  in  accuracy  of  simplified  transformations  between  natural 
and  physical  coordinate  systems  may  become  unacceptable.  With  stresses  assumed  in  physical  coordinates, 
the  present  study  does  not  consider  minimum  expansions  of  assumed  stress  modes  and  the  expansions  are 
assumed  complete  to  the  highest  order  present  in  the  strain  field.  Although  this  will,  in  cases,  result  in 
significantly  greater  number  of  stress  terms  than  the  minimum  required  to  suppress  zero  energy  modes,  the 
completeness  property  of  the  assumed  stresses  preserves  invariance  and  avoids  spurious  kinematic  modes  as 
the  vector  space  formed  by  the  assumed  stresses  is  guaranteed  to  span  the  strain  space.  In  the  use  of  the 
Hellinger-Reissner  functional,  stress  constraints  need  not  be  enforced  a  priori,  however,  they  can  be  applied  to 
satisfy  field  equilibrium  and  compatibility  conditions  pointwise  in  order  to  reduce  the  number  of  independent 
stress  modes.  In  physical  coordinates,  element  invariance  will  be  preserved  under  field  operations  of  elasticity 
if  complete  expansions  for  the  stresses  are  used.  Stress  expansions  assumed  in  natural  coordinates  permit 
greater  freedom  in  the  selection  of  expansions  for  specific  stress  components.  This  flexibility  allows  a  degree  of 
tailoring  of  element  strain  energy  mode  representation  while  maintaining  invariance  and  a  reduced  sensitivity 
to  mesh  distortion. 


5  Determination  of  Explicit  Forms 

A  procedure  for  simplifying  the  expressions  involved  in  (11)  by  utilizing  permissible  transformations  of 
assumed  stress  fields  is  now  detailed.  A  two-step  transformation  of  the  assumed  stresses  is  suggested  by  an 
examination  of  the  flexibility  matrix  given  in  equation  (6)  written  out  fully  as 

H  =  f  [\J\pT’SP]d(,dridC.  (13) 

An  initial  observation  is  that  the  structure  of  the  integrand  in  (13)  can  be  simplified  through  an  apportioning 
transformation  of  the  material  compliance  matrix,  5,  to  the  stress  modes  in  P  thereby  allowing  further 
simplifications  in  the  flexibility  matrix  to  be  achieved.  Towards  this  aim,  the  assumed  stresses  are  first 
transformed  through  the  introduction  of  a  symmetric  ’distributing’  matrix,  D,  which  acts  to  subsume  the  5 
matrix  into  P  via  an  identity  operator  as 

P  =  IP  =  (DD-^)P  =  D(D-^P)  =  DP  (14) 

where  the  distributing  matrix  is  defined  as 

D  =  5-*/-  (15) 

Although  formally  permissible,  this  operation  transforms  the  initially  assumed  stress  modes  into  a  set  of 
vector  polynomials,  P,  which  do  not  have  a  direct  physical  interpretation.  Instead  of  introducing  new 
terminology,  however,  the  transformed  stresses  in  P  will  simply  be  referred  to  as  stress  modes.  The  inverse 
squaure  root  of  the  compliance  matrix  is  obtained  via  a  standard  spectral  decomposition  as 

S  =  QAQ^  (16) 

in  which  Qis  a  column  matrix  of  the  normalized  eigenvectors,  of  S  given  by 

which  is  a  unitary  matrix  with  the  property  that 

Q^  =  Q-' 
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and  A  18  a  diagonal  matrix  formed  by  the  eigenvalues,  <pi,  of  S  given  by 

A  =  diag[(fii,ip2,  —,'Pn] 

where  n  =  dim{S).  With  the  above  definitions,  the  D  matrix  is  given  by 


D=  5"*/*  =  QA-*/’(3^  (17) 

where 

The  symmetry  and  positive  definiteness  of  the  material  property  matrix  guarantees  the  existence  of  the 
decomposition  and  explicit  expressions  for  Q  and  A  for  both  2-D  and  3-D  orthotropic  compliance  matrices 
are  presented  in  Appendix  I. 


Substitution  of  (14)  into  (13)  yields 

7/  =  /  J  J  [\J\P'^D'^SDP]d^dtjdC 

where,  from  the  definition  of  D  and  the  symmetry  of  both  S  and  /?,  we  obtain 

D^SD  =  =  SS"*  =  I 

and  the  flexibility  matrix  reduces  to 

ri 


//=  j  j  J  [\J\P'^ PHdrtdC 


(18) 


(19) 


(20) 


A  second  field  transformation  is  motivated  by  the  form  of  equation  (20)  which  suggests  its  use  to  define  an 
inner  product  space  where  a  Gram-Schmidt  procedure  can  be  employed  to  generate  an  orthonormal  spanning 
set  of  stress  modes,  P*,  which  are  a  special  linear  combination  of  the  modes  present  in  P.  The  weighted 
inner  product  is  therefore  defined  as 


<  Pi,Pj  >=  y*  PAd^drjdC  =  6ij 


(21) 


where  6ij  is  the  Kronecker  delta  function.  The  linear  combination  yielding  a  sequence  of  orthogonal  stress 
modes  is  defined  by 


Vi  =  Pi-'^<  P’,Pi>  pr 

}=i 

which  are  normalized  to  form  basis  vectors,  P/,  as 

pr  =<  Vi,  Vi  v; 


(22) 


(23) 


Substitution  of  P*  into  equation  (20)  yields  by  definition 

H  =  j'j'j'[\J\P-'^P-]d^drfdC  =  I  (24) 

Hence,  by  fully  exploiting  permissible  operations  on  the  assumed  stress  modes,  the  transformations  due  to 

the  ^plication  of  a  distributing  matrix  and  the  generation  of  a  weighted  orthonormalized  basis  yield  the 
element  stress  field  as 

P  =  Z?P*  (25) 

and  the  flexibility  or  complementary  energy  matrix  ’H’  is  eliminated  by  formally  reducing  it  to  a  matrix 
identity.  The  expression  for  the  element  stiffness  matrix  is  now  given  by 

I<  =  G^G  (26) 
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wh«e 


(27) 


G  =  [I  -  R{HT R)-^ R'^]G 

Separating  out  the  Jacobian  determinant  from  the  compatible  and  incompatible  strains  as 


B  =  (iA^)  =  i7|fl*(^.»?.0 

(28) 

fl  =  (LM).=  |lfl-(€,»,,0 

(29) 

and  substituting  into  (7)  and  (8),  the  constituent  matrices  become 

(30) 

R=  j'j'j'y'^DB']didr)dC 

(31) 

With  the  removal  of  the  flexibility  matrix,  the  element  stiffness  matrix  is  fully  determined  by  the  two 
constituent  matrices  in  (30)  and  (31)  which  represent  the  elastic  strain  energy  contributions  of  the  assumed 
stresses  and  the  compatible  and  incompatible  strain  modes.  The  absence  of  the  Jacobian  determinant  in 
the  denominator  allows  a  direct  computation  of  linear  algebraic  forms  for  the  G  and  A  matrices  based  on 
the  regular  structure  of  the  strain  modes  and  the  transformed  stress  fields.  Explicit  expressions  for  various 
element  matrices  will  be  developed  in  subsequent  sections.  In  the  computation  and  manipulation  of  element 
matrices,  however,  not  all  operations  are  most  efficiently  obtained  explicitly  and  numerical  procedures  will  be 
prescribed  for  certain  procedures  when  computationally  advantageous.  For  example,  equation  (27)  requires 
the  inversion  of  an  inner  matrix  product  in  the  definition  of  G  given  by 

A  =  Z-' {R^ R)-'  (32) 

However,  because  the  columns  in  A  are  equal  to  the  number  of  incompatible  modes,  the  dimension  of  A  is 
given  by 

dim(A)  =  dim(X)  x  dim(A)  (33) 

which  is  usually  small.  In  addition,  because  A  is  symmetric,  the  inversion  effectively  reduces  to  the  com¬ 
putational  effort  of  inverting  a  matrix  of  triangular  form  which  can  be  computed  explicitly  in  genered  but, 
with  increasing  matrix  order,  a  numerical  scheme  is  preferred.  A  second  example  is  the  Gram-Schmidt 
procedure  which  quickly  leads  to  cumbersome  expressions  for  the  coefficients  of  the  orthogonal  stress  modes 
when  performed  symbolically.  However,  with  an  initial  evaluation  of  the  basic  scal»  integrals  required  in  the 
weighted  inner  product  involving  the  Jacobian  determinant  and  powers  of  the  assumed  stress  polynomials,  a 
simple  numerical  procedure  may  be  used  for  efficiently  computing  the  linear  combination  of  modes  present 
in  P  to  generate  the  orthonormalized  basis  set  P*.  Symbolic  representations  will,  however,  be  generated  for 
selected  elements  used  in  the  present  study. 

The  explicit  form  of  the  element  stiffness  matrix,  decomposed  into  contributions  due  to  compatible  and 
incompatible  displacements,  is  given  by 

~  ^ in  ~  9ni9nj  ~  Jtj  '  ^}i  —  ^ij  (34) 

where  the  components  of  G  and  A,  gtj  and  Vij,  are  obtained  by  integrating  (30)  and  (31)  and  Oij  are  the 
components  of  the  A  matrix  given  in  equation  (32).  The  various  indices  range  from 


t  =  1,2,3..., fV,;  j  =  1,2,3,. ..,i;  n.fc  =  1,2,3..., AT^;.  s,m  =  1,2,3..., Na 


where  denotes  the  number  of  element  degrees  of  freedom,  Np  number  of  independent  assumed  stress 
nnodes,  and  Nx  the  number  of  incompatible  displacement  modes. 

The  above  method  for  determining  explicit  forms  of  element  stiffness  matrices  will  be  demonstrated  with 
4-node  quadrilateral  plane  and  $-node  hexahedral  solid  element  formulations.  The  explicit  integration  of 
(26)  offers  a  significant  decrease  in  computational  cost  over  a  purely  numerical  evaluation  of  (11). 
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6  4-Node  Plane  Quadrilateral  Elements 


Figure  1.  Quadrii&teral  element. 


The  stiffness  matrices  of  several  different  4-node  plane  elements  will  be  explicitly  derived  in  this  section. 
Elentent  configuration  and  node  numbering  is  shown  in  Figure  1 .  Two  complatable  elements  are  presented 
to  highlight  features  of  both  the  developed  procedure  and  the  hybrid  stress  method  followed  by  two  ele¬ 
ment  formulations  incorporating  incompatible  displacement  modes  to  optimize  element  performance.  Stress 
expansions  are  assumed  in  both  physical  and  natural  coordinates  to  denK>nstrate  the  generality  of  the  de¬ 
veloped  methodology.  The  correspondence  to  existing  element  formulations  will  be  identified. 

The  displacement  functions  u,  are  given  by 


The  mapping  from  local  physical  coordinates  to  natural  coordinates  is  given  by 

■c  =  ai4  +  «2»;  + 
y  =  63fi7 

where 


Ufl  ho 

■  1  1  1  1  ■ 

yi 

ai  hi 

_  1 

-1  1-11 

X2  yz 

<l2  bn 

“  4 

-1-1  11 

*3  ys 

O3  63 

1-1-11 

.  *4  y4  . 

(35) 


(36) 


(37) 


The  compatible  Element  E4PQ 


A  4-node  plane  element,  designated  E4PQ,  is  explicitly  derived  using  stresses  assumed  in  physical  coordi¬ 
nates.  The  stress  field  is  selected  as  complete  linear  expansions  with  the  constrzunt  <t  =  0  applied  resulting 
in  the  following  seven  equilibrating  stress  modes 


P  = 


The  derivation  of  weighted  orthonormalized  stress  modes  are  obtained  in  the  form 


Pi  I  +  Pi2X  +  P, 

PiA  +  PiiX  +  piey 
Pii  +  Pis-c  +  Pioy 


i3y  1 
•  6J/  > 

•9y  J 


and  are  given  by 


P"  = 


Pii 

0 

0 


0 

P2A 

0 


0 

0 

P37 


P^l  +  P43j/ 


P51  +  P52*  +  P53y 
0 


P61  +  P62*  +  P63y 
P64  +  Pssy 


P7i  +  P72X  +  pray 
II  d  P64  +  Pssy  P7A  +  P7iX  Prsy 

P47  +  PAaX  P57  +  P5SX  psgy  p67  +  PSS*  +  PSSy  P77  +  PTSX  +  P79y 


(38) 


(39) 


(40) 
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The  coefficients  of  the  stress  modes,  pij,  are  presented  in  Appendix  II.  Because  the  stress  modes  in  (38)  are 
self  equilibrating,  equilibrium  is  unaffected  through  the  linear  combination  of  modes  leading  to  (40). 

The  integration  of  equation  (26)  is  obtained  in  a  straightforward  fashion.  Using  the  following  constants 


«!* 

N 

CO 

II 

-*1^2 

64* 

= 

—  032* 

62* 

=  *3^1 

-  f>iZ^ 

6S* 

=  ai^3 

-0321 

where 

63* 

=  *2^3 

-  ^3^2 

cst 

=  032* 

-022! 

the  components  g„t  of  (30)  are  given  by 


T 

^3 

1 

-1 

-i 

1 

2 

1 

-1 

-1 

3 

-1 

1 

-1 

4 

1 

1 

1 

ffn(2t-l)  =  Cn.(dnPnl  +  djoPna)  +  d33p„7C4|;  +  ^{(Pn’^fll  +  Pn5di2)(aie2i  +  02C3t)+ 
(Pnsdll  +  Pn6dl2){6ie2t  +  +  d33P„8(aiC5t  +  <J2e6t)  + 

<f33Pn9(ilC5t  + 


ffn(2t)  =  e4k(dnp„i  +  d-ilPnA)  +  ds^Pn'^lk  +  5{(Pn2dl2  +  Pnlid22){<^\^Sk  +  0266*)  + 
(Pn3dl2  +  Pn6d22)(6lC51'  +  f>2«6t)  +  d33P„8(aie2t  +  02^3*)  + 
d33Prt9(6ie2t  +  b-tezi:)} 


(41) 


(42) 


where  dij  are  elements  of  the  distributing  matrix.  The  stiffness  matrix  is  then  given  explicitly  as 


—  9’i9nj 


i=l,2,3 . 8;  j=  1,2,3,  ...,j; 


n=  1,2,3,  ...,7 


(43) 


and 


Kji  =  Kii 


The  stiffness  matrix  given  above  is  equivalent  to  the  7-/?  hybrid  element  using  the  assumed  field  in  (38) 
and  is  closely  related  to  the  standard  minimum  5-/?  hybrid  element  of  Pian  which  incorporates  equilibrating 
stress  expansions  given  by 


P  = 


1  0  0  y  0 

0  1  0  0  X 

0  0  10  0 


(44) 


A  validation  of  the  procedure  is  presented  in  Table  1  by  showing  the  equivalence  of  eigenvalues  obtained 
from  the  explicit  element  stiffness  matrix  given  by  (43)  and  the  numerical  computation  using  equation  (11). 
For  generating  the  results  in  Table  I,  a  general  quadrilateral  configuration  was  arbitrarily  selected  as  shown 
in  Figure  2.  A  2  x  2  Gaussian  integration  rule  was  used  for  the  numerical  stiffness  matrix. 


(0^  1.0) 


Fiqure  2.  General  quadrilateral  element. 


Table  1.  Equivedence  of  eigenvalues  obtained  from 
the  explicit  and  numerically  integrated 
form  of  the  element  stiffness  matrix. 


A 

Numerical  K  matrix 

Explicit  K  matrix 

1 

0.0 

0.0 

2 

0.0 

0.0 

3 

0.0 

0.0 

4 

500.69450 

500.69450 

5 

585.41852 

585.41852 

6 

1190.1325 

1190.1325 

7 

1297.3502 

1297.3502 

8 

2172.8470 

2172.8470 
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The  compatible  Element  E4PR 


The  E4PR  element  is  formulated  using  unconstrained  stress  expansions  resulting  in  the  following  nine  inde¬ 
pendent  stress  modes 

[izyOOOOOO" 

P=  OOOliyOOO  (45) 

000000  IxyJ 

In  performing  the  initial  apportioning  transformation  of  the  assumed  stress  heiii  as  given  by 

P=D-‘P  (46) 

the  stress  modes  become  coupled  due  to  the  action  of  the  distributing  matrix  which  complicates  the  subse¬ 

quent  orthonormalisation  procedure.  However,  a  simplification  can  be  made  through  a  linear  combination 
of  the  modes  in  (46)  with  constant  terms  being  absorbed  into  the  vector  of  unknown  expansion  coefficients, 
0.  With  stresses  expressed  by 

a  =  Pp=  DP0  (47) 

the  modes  may  be  rearranged  and  the  P  matrix  decomposed  to  give 


1  0  0 


0  0  1 


D-*<t  =  D-‘  0  1  0  l3i  +  D-^\0  X  0  /32  +  D- 


X  0  0 


0  0  X 


y  0  0  ■ 

0  y  0  ^3 

0  0  y 


where  0i  are  subvectors  of  0.  Rewriting  the  above  as 

D-^<t  =  D-^I0i+  xD-^I0i  -t-  yD-‘//?3 
suggests  a  linear  combination  of  the  stress  modes  for  each  partition  defined  by 

0i  =  D-^0i 

which  leads  to  the  simplification 

D~‘<r  =  D~^ I D0i  +  xD~^ ID02  +  yD~^  ID03  ~  I0i  +  xI02  +  yI03 


and  equation  (47)  becomes 


<T  =  DP0 


Although  for  arbitrary  stress  expansions  equation  (14)  is  strictly  valid,  because  the  stress  expansions  defined 
in  (45)  are  balanced  for  all  components,  the  inverse  distributing  matrix,  D~^,  can  be  completely  removed 
and  the  form  of  P  is  made  identical  to  P.  With  0  defined  by  equations  (9)  and  (10),  the  linear  operations  on 
P  are  automatically  accounted  for  and  the  distinction  between  0  and  0  may  be  neglected.  Orthogonalizing 
the  assumed  stress  is  thus  reduced  to  determining  an  orthonormal  sequence  of  scalar  functions  using  the 
weighted  inner  product  defined  in  (21)  and  the  basis  functions  given  by  [I,x,y]. 

The  weighted  orthonormalized  stress  modes  are  obtained  as 

"  pI  Pj  P3  0  0  0  0  0  O' 

p*  =  0  o'  0  p;  p5  p5  0  0  0  (53) 

0  0  0  0  0*  0  pj  p5  p5 


where 


Pi  =  Pii 

PZ  =  P2I  +  P22* 

P3  =  P3l  +  PZlX  +  P33y 


The  coefficients,  pij,  are  presented  in  Appendix  II. 

The  integration  of  (30)  yields  the  components  gtj  with  i  =  1, 2, 3  smd  j  =  1, 2, 3, 4  as 


ff(iK2i-i)  = 


<7(i+3)(2;-l)  =  di2^i; 

?(*+3)(2i)  = 


ff(«+6)(2i-l)  —  <^33i/’ij 
ff(i+6)(2j)  =  d33<l>ij 


where 


’+  [pi2(oieji  +  a^esj)  +  pisibicy  +  62e3j)]/3 
V’ij  =  Pii«4j  +  (piaC^i^si  +  a3«6j)  +  PisibiCsj  +  ^e6j)j/3 

and  where  dij  are  elements  of  the  distributing  matrix  and  e,-;  are  defined  in  (41). 
The  stiffness  matrix  is  explicitly  given  by 


^ij  —  Smiffmj 


i  =  1,2,3....,8:  j  =  1,2,3, 


i;  m=  1,2, 3,. ..,9 


The  interest  in  this  hybrid  formulation  lies  in  the  limitation  principle  of  [11]  which  states  that,  in  the  limit 
of  the  order  of  unconstrained  stress  expansions,  the  hybrid  stress  method  converges  to  the  stiffness  matrix 
obtained  from  a  purely  displacement-based  formulation.  For  a  parallelogram,  the  hybrid  E4PR  element 
exactly  duplicates  the  displacement-based  stiffness  matrix.  It  is  of  interest  to  quantify  how  well  the  explicit 
hybrid  element  formulation  can  be  substituted  for  the  displacement-based  element  which  does  not  permit 
a  simple  integrated  representation  under  general  element  distortion.  Figure  3  shows  an  initial  unit  element 
geometry  and  its  change  as  a  funtion  of  a  distortion  parameter,  e.  Table  2  presents  ratios  of  traces  of  the 
element  stiffness  matrices  of  increasingly  distorted  element  configurations  as  a  function  of  the  distortion  pa¬ 
rameter  comparing  the  hybrid  element  with  the  corresponding  displacement-based  element.  A  2  x  2  gaussian 
quadrat.'ire  rule  was  used  in  computing  the  displacement  -  based  matrix.  It  is  shown  that  the  explicit  hy¬ 
brid  o.re88  formulation  yields  a  consistently  more  flexible  element  and  offers  a  clear  computational  advantage. 


Ratios  of  stiffness  matrix  traces  comparing  explicit 
hybrid  method,  K/t,  with  numerical  displacement-based 
method,  Kj  under  increasing  element  distortion. 


e 

Trace(Kx)/Trace(Ki) 

0.0 

1.000000 

0.1 

0.998322 

0.2 

0.993153 

0.3 

0.984084 

0.4 

0.970401 

0.5 

0.951032 

0.6 

0.924311 

The  Incompatible  Element  E^PL 


The  E4PL  element  incorporates  the  unconstrained  stress  modes  given  by  equation  (45)  with  the  addition  of 
incompatible  displacement  modes.  The  incompatible  displacement  functions  are  selected  as  17^]  which  are 
then  modified  according  to  Reference  12  to  identically  satisfy  the  strong  form  of  the  convergence  requirement 
on  incompatible  modes  given  by 

r 

Luxdv  =  0 


I 


which  yields 


r«i 

1  _ 

1 « J 

Ml 

0 

Mo 

0 

A: 

Aj 

0 

Ml 

0 

M2 

< 

As 

.  . 

where 

Mi  =  C--UM  +  f2r}) 

Ml  =  -I-  |(/i^  +  fivi) 

By  integrating  (31),  the  elements  of  the  R  matrix  with  n  =  1,2,3  are  given  by 


(57) 


(58) 


(59) 


rn.i  =  diie‘ 

»‘n.2 

=  dioBl 

»'n,3 

=  du0n 

»'n,4  =  di20® 

^n+Z.i  =  di20i 

'■n-^3,2 

=  ^2202 

^n+3.3 

=  d^ex 

’’n-l-3,4  =  ^220® 

(60) 

=  dssSj 

'■n-(-6,2 

=  d330^ 

rn+6,3 

=  d330^ 

Tn-fS.^  =  ^330® 

where 


Qn  =  +  Pn2h2  +  Pn3A3)/3  0*  -  +  Pnll^a  -  Pnal^e)/^ 

Qn  =  4(pnlh4  +  P„2A8  +  Pn3h6)/3  6®  =  4(p„l/l9  +  P„2hl0  —  Pn3/>6)/3  (61) 

©n  =  4(p„i/»7  +  p„2h3  +  Pn3h3)/3  9®  =  i(Pnlf^7  +  Pni^^  -  PnsAsj/S 

where  the  constant  terms  ft  and  hi  are  presented  in  Appendix  III.  The  elements  of  the  Z  matrix  defined  in 
(32)  are  computed  as 

—  fni^nj  i  2ji  —  Zij  (fi2) 


where 


1  =  1, 2, 3, 4;  i  =  l,2, 


i;  n  =  l,2,3 . 9 


The  components  a,j  of  the  inner  product  defined  in  (32)  are  obtained  through  a  closed-form  inversion  of  the 
symmetric  A  x  4  Z  matrix.  The  element  stiffness  matrix,  A',  is  given  by  the  sum  of  contributions  due  to 
compatible  euid  incompatible  displacements 


where 


ftij  —  A^,j  +  —  9ni9nj  ~  9ni^n$f^$m^km9k}  i  hCji  —  Kij 

i=  1,2,3,...,8  ;  j  =  1,2,3,..., i  ;  n,fc  =  1,2,3,...,9  ;  s,m=l,2,3,4 


The  Incompatible  Element  E^NL 


The  E4NL  element  is  formulated  using  unconstrained  stress  expansions  assumed  in  natural  coordinates 
resulting  in  the  following  nine  independent  stress  modes 


] 

■  1 

n 

0 

0 

0 

0 

0 

0  ■ 

= 

0 

0 

0 

1 

»? 

0 

0 

0 

/ 

0 

0 

0 

0 

0 

0 

1 

7 . 

01 

02 


)=T0 


09  ) 


(63) 


In  order  to  preserve  the  constant  stress  terms,  the  natural  stresses  are  mapped  to  physical  space  through  a 
contravariant  transformation  using  Jacobians  computed  at  the  element  centroid 


or 

The  initial  stresses  are  transformed  as 
where 

and  the  constants,  Oj,  are  given  by 

=  diiai-fdi26i 

0t2  —  012^1  +  022^1 

Qt3  =  <1330102 


=  (4)?(4)<r‘^ 


p  =  nr 


p  = 

D-^T„r 

=  fr 

ai 

04 

07 

t  = 

02 

05 

<»8 

03 

06 

O9 

04 

=  <liio§  -b  <11262 

07 

as 

^  ^13^2  ^22^2 

Os 

as 

=  <^3361 62 

O9 

(64) 

(65) 

(66) 

(67) 

2(<fii0i02  +  0l2hll>2) 

2(dl2®l®2  "b  022h\l>2)  (68) 

<133(0162  -b  0261) 


in  which  0ii  are  elements  of  the  inverse  distributing  matrix, 
modes  defined  by 

r  f 


0= 


T 


T 


By  introducing  a  linear  combination  of  stress 
0  (69) 
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the  coupling  of  the  transformed  stresses  can  be  eliminated  resulting  in 

’l^ijOOOOOO' 

P=  0001^i;000  (70) 

OOOOOO’l^i? 

and  the  orthonormalization  is  effectively  reduced  to  determining  an  orthonormal  sequence  of  scalar  functions 
using  the  basis  set 

The  weighted  orthonormalized  stress  modes  are  obtained  as 

■  pI  p5  p3  0  0  0  0  0  O' 

p*  =  0  0  0  pI  p;  0  0  0  (7i) 

0  0  0  0  0*  0  Pi  p;  p5 


where 

Pi  =  Pii 

Pj  =  P21  +  P22f 

P3  =  P31  +  P32^  +  P33»7 

and  the  coefficients,  p,j,  are  presented  in  Appendix  III. 

The  integration  of  (30)  yields  the  components  gmn  with  m  =  1,2,3  and  n  =  1,2,3, 4  as 

jf(m)(2n-l)  =  du^mn  jf(m+3)(2n-l)  =  dl2^mf»  P(m+6)(2n- 1)  = 

i7(m)(2n)  =  di2V>mn  ?(fn+3)(2n)  =  9(m+6)(2n)  = 

where 

0mn  •  PmlCln  +  \Pm2^in  + 

V’mn  —  Pml^An  +  (Pm2^5n  +  Pm3*6n]/3 


(72) 

(73) 

(74) 


and  where  dij  are  elements  of  the  distributing  matrix  and  e,'i  are  defined  in  (41). 

The  incompatible  displacement  modes  are  the  same  as  those  used  in  E4PL  and  yield  the  components  r^n 
of  (31)  with  m—  1,2,3  as 


^m,l 

=  diiei. 

rm.7  =  dijej. 

rm.3  =  dll©;!, 

^m,A  =  di20m 

^m+3,1 

=  di2©i. 

rm+3,2  =  d22©m 

rm+a.S  =  di2©m 

»‘m+3,4  =  d220m 

(75) 

*m+6,l 

=  d33©5. 

Tm+e.r  =  ^330^ 

»‘m+6,3  =  d330^ 

»’m+6,4  =  d330^ 

where 

©m  =  Pmlhl  +Pm2h2  +Pm3h3  ©m  =  Ptnl^J  +  Pmzf^a  —  PmzfiZ 

©m  =  PmlhA  +  Pmlhi  +  Pm3f*6  ©m  =  Pml *9  +  Pm2hl0  -  Pmsh-r  (76) 

©m  =  Pmlh? +  Pm2/>2 +Pm3h3  ©m  =  Pmlhj  +  Pmzf^a  -  Pmzf^l 

The  constant  terms  ft  and  hi  are  contained  in  Appendix  III.  The  elements  of  the  symmetric  Z  matrix  are 
computed  as 

where 

1=1,2, 3, 4;  j  =  l,2,...,i;  n  =  1,2,3, ...,9 


The  components  atj  are  then  obtained  through  the  inversion  of  the  symmetric  4  x  4  Z  matrix.  The  element 
stiffness  matrix,  K,  is  given  by 

—  ^^q,j  4"  —  9ni9nj  ~  9ni^nt^$m'^km9kj  i  —  diij 


where 


i=  1,2,3,...,8  ;  i  =  1,2,3,..., »■  ;  n.ifc  =  1,2,3,...,9  ;  s,m=l, 2,3,4 
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The  E4PL  and  E4NL  elements  are  related  to  the  Pian-Sunihara  element  [13]  which  has  demonstrated  excep- 
tionrd  performance  in  plane  stress/plane  strain  problems.  An  elegant  derivation  of  an  explicit  form  for  the 
Pian-Sumihara  element  is  presented  in  Reference  [14]  utilizing  a  scaling  procedure  and  stabilization  matrices 
to  obtain  stiffness  coefficients.  The  present  method  avoids  the  need  for  numerical  stabilization  and  provides 
a  generic  approach  for  obtaining  explicit  element  stiffness  matrices  in  hybrid  element  formulations. 


7  8-Node  Hexahedral  Continuum  Elements 


Figure  4.  Hexahedral  element  configuration. 

Two  8-node  hexahedral  continuum  elements  incorporating  incompatible  displacements  will  be  considered 
in  this  section.  Element  geometry  and  node  numbering  convention  is  shown  in  Figure  4.  The  first  element, 
designated  E8PL,  is  based  on  unconstrained  stress  expansions  assumed  in  physical  coordinates.  The  second 
element,  designated  E8NL,  is  based  on  stresses  assumed  in  natural  coordinates.  The  selection  of  incompati¬ 
ble  modes  are  identical  to  those  used  in  [12], 


The  displacement  functions  Uf  are  given  by 

“*  =  I  I  =^^(l+^<^)(l  +  »;.»7)(l+CiC)|  Vi  |=fV<9 

The  isopsirametric  mapping  between  local  physical  and  naturetl  coordinates  is  given  by 

X  =  ui^  -f  a^t]  +  uaC  -f  a4^r)  +  +  06»7C  +  07^»/C 

y  =  +  hi  +  ^sC  +  +  hlC  + 

z  =  ci^  4-  ciT]  +  caC  +  c^jj  +  Cs^C  +  c^lC  +  crklC 

where 
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1  -1  -1 

1-1  11 

-1 

Xo 

ye 

76 

06 

bo 

C6 

1  1  -1  - 

-1  -1  -1  1 

1 

X7 

yr 

77 

07 

67 

C7 

-1  1  -1 

1  1-11 

-1 

. 

ys 

^8  . 
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(77) 


(78) 


(79) 


The  E8PL  Element 


For  an  S-node  isoparametric  continuum  element  with  trilinear  displacement  interpolants,  the  number  of 
independent  strain  modes  is  18.  A  detailed  study  of  strain  modes  in  Reference  [9]  has  determined  that  a 
linear  expiuision  field  is  insufficient  and  will  give  rise  to  spurious  zero  energy  modes  unless  sefected  quadratic 
terms  are  added.  In  physical  coordinates,  using  the  completeness  property  to  guarantee  element  invariance, 
the  number  of  assumed  str^  modes  is  60.  Equilibrium  and  compatibility  constraints  can  be  appled  to  reduce 
the  independent  modes  to  42  but  is  not  performed  in  the  present  analysis.  Several  versions  of  the  E8PL 
element  are  formulated  using  different  assumed  stress  fields  to  assess  the  effect  on  element  performance.  In 
E8PLi,  stresses  are  assumed  as  complete  linear  expansions  for  all  stress  components.  This  selection  yields 
an  element  containing  2  zero  energy  modes  which  would  have  to  be  removed  through  stabilization  to  yield  a 
useful  element.  Stress  expansions  including  quadratic  cross  terms  are  incorporated  into  the  E8PL2  element. 
This  element  possesses  the  requisite  number  of  rigid  body  modes  but  is  not  invariant  under  coordinate  trans¬ 
formation.  A  third  element,  E8PL3,  incorporates  complete  quadratic  expansions  which  precludes  spurious 
kinematic  modes  and  ensures  invariance.  This  element,  however,  is  shown  to  suffer  from  significant  sensi¬ 
tivity  to  distortion.  The  three  versions  of  the  E8PL  element  are  compared  to  a  3-D  element  formulated  in 
natural  coordinates  which  demonstrates  optimal  behavior. 

To  encompass  the  various  stress  fields  assumed  in  the  E8PL  element,  the  initial  stress  field  is  selected  as 


P  = 


[A] 


(A1 


[Pi] 


[Pi] 


[Pi] 


[Pi]\ 


which  are  complete  quadratic  expansions  given  by 

[Pi]  =  [l,I,y,^,a:y,y^,rx,x^y^r^] 


Orthonormalizing  yields 

[Pi]  =  [Pi.P2iP3>P4>P5>P6>P7<P8>Po>PIo] 
where,  for  t  =  1,2,3,...,  10,  the  general  form  of  each  mode  is  given  by 


(80) 


(81) 

(82) 


P< 


Pii 


+  Pi2*  +  p-ay  +  PhZ  -I-  pJsiy  +  p^yz  -|-  pljzx  -I-  H-  pjgy^  + 


(83) 


and 

Pit  =  0  for  k  >  i 

A  procedure  for  generating  the  constant  coefficients  in  the  expressions  for  p,*  is  presented  in  Appendix  II. 


The  stress  fields  are  selected  in  the  following  expressions  by  setting  Ng  equal  to  4,7  and  10  for  the  E8PL1, 
E8PL3  and  E8PL3  elements  respectively.  Integrating  equation  (30),  the  components  gmn  where  m  = 
1,2, 3, ...,  TVg  and  n  =  1,2, 3,...,  8  result  in  lengthy  linear  eilgebraic  expressions  which  are,  however,  highly 
structured  and  allow  a  compact  notation  to  be  used.  The  expressions  for  gmn  are  given  by 


ff(mK3n-2) 

flf(m)(3n-l) 

9(m)(3n) 

0(m+iV.)(3n-2) 

ff(m+7V.)(3n-l) 

9(m+N,)(an) 


= 

dliei^n 

P(m+2/V.)(3n-2) 

= 

‘^310mn 

)(3n  -  2) 

= 

= 

ff(m+2//,)(3n-l) 

d320mn 

9(m+AN,)(an-l) 

= 

0 

= 

dlzQ^mn 

9(m+2N.)(3n) 

=  _ 

d33ef,„ 

9(m+AN.)(3n) 

= 

‘^550mn 

= 

dn^mn 

^7(m+3^.)(3n-.2) 

0 

P(m+5JV.)(3n-2) 

deeOl^ 

= 

^220^ 

ff(m+3yv.)(3n-l) 

dA^e^mn 

9(m+iN.)(3n-l) 

d660mn 

= 

*^230mn 

ff(m+3Ar.)(3n) 

= 

ff(m+5JV,)(3n) 

= 

0 

(84) 
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where 


eL  =  J^p;„,^Jn  (85) 

«sl 

For  m  =  1 

0t,  =  -wii/3  (86) 

For  m  =  2,3,4 

^ni  =  -(3ait«a,  +  3aju^  +  +  a"ti>5i  +  aj  u;*,  +  aju>7i  +  2a^w^)l21  (87) 

For  m  =  5, 6, 7, 8, 9, 10 

=  -[«?( 10/3? «’«  +  /??  *"9.  +/37«'ia+/3ru'i3i+/34“'ioi  +  5/^»"6i  +  5/3Ju;J,+ 
o5(10/^a»^  +  0lw\ai  +  +  0^X0^^^+ 

O5(10/?y  tl>|,.  +  +  06^21i  +  5/?f  +  J®?  “»l«  + 

oj(10^u»|i  +  +  0S^i2i  +  /3yu>f«  +  /3i  u^oi  +  /3?  (88) 

«8  (lO^tuij  +  /}^w^  +  /3?u;fj,  +  +  /3?«'j2i  +  /??  «»L  +  05^36$+ 

a2(lO0?wl  +  /??«>{„  +  06»>iai  +  /^*"2U  +  04'^23i  +  /35  «"24.  +  /^T  «'27i+ 

^"(/^y  +  07^l9i  +  05^22i  0?^^  +  03^244  +  ^*‘'j7i]/185 


The  constants  Wj^  are  given  by 

Wii  =  (is4  +  34jj)*i  +  (fiea  +  ifs)*j  +  (^m  +  i*i  )^3  +  (^Je  +  +  (ijj  +  f34)se  +  {6^^  +  3<ji  )*) 

Wji  =  3(634  +  653)*!  +  (672  +  3613)23  +  (637  +  3631)23  +  (653  +  634)24  +  (647  +  3613)25  +  (675  +  3641)27 

WSl  =  (^74  +  3633)21  +  (3633  +  6^7)23  4-  (633  +  6*1  )24  +  3(643  +  ^!«)*l  +  (^73  +  36*1)2*  +  (6*7  +  3634)27 

w*!  =  (6^7  +  36m)2i  +  (36*3  +  ^7j )*3  +  (^62  +  ^fs)*4  +  (^7*  +  36^)2^  4-  (36*3  +  ^*7)*«  +  3(6^1  4-  635)2) 

xx>ii  =  3(^n  +  +  3(6m  4-  6^5)2)  4-  (6^7  4-  36^4)2)  4-  (6^5  4-  36*1)2)  4-  3(6y5*4-  6^7)2)  4-  (6*5  4-  36^  )2# 

«*»•  =  3(^37  +  6m)«)  +  (67*  4-  36^3)23  4-  3(6*4  +  ^ei)*3  +  (36^3  4-  6^4)2)  4-  (6^7  4-  36f5)2)  4-  3(6^5  4-  6*1  )27 

wf«  =  (^$7  +  365j)23  4-  (6*4  4-  36^2)23  4-  (36*3  +  ^4s)*4  +  3(6^*  4*  6*3)2)  4»  3(6^  4-  6^)2)  4-  3(6**  4-  637)2) 

■tn*j  =  6m^  +  6^*2)  4-  6^52) 

40*i  =  (9657  4-  1563*)2j  4-  (96*1  +  156*3)23  4-  (56*2  4-  8655)2)  4-  5(6**  4-  3653)2*  4-  6(637  4-  8653)2*  4-  15(65i  +  6*5)2) 

Wiot  —  5(657  +  363* )2i  4"  5(671  +  36*3)23  4-  (96*3  +  5655)24  4"  (96**  4"  15653)2*  4"  (86*7  4-  1565s)2*  4"  15(6*i  4-  635)2) 

Will  —  (3657  +  563* )2i  4-  (36*1  4-  56*3)23  4-  3(6*3  4-  65s)24  +  (36**  4-  56^)2*  4-  (8637  4-  5653)2)  4-  5(6*i  4-  635)2) 

Wi3j  =  (36*7  4"  5635)23  4"  (36*4  +  56*3)23  4-  (36**  4"  56*3)24  4-  5(6^  4-  6*3)2*  4"  5(6*3  +  55s)2*  4-  5(6*4  +  5*7)*)  ( 

=  3[(9654  4-  15652)*;  +.  (56**3  +  965*)*)  4-  (565*  4-  964*1  )*)  4  5(6)*  4-  3653)*)  +  5(65*3  +  «54)*i  +  5(6***  4-  363*,  )*)] 
w54i  =  3[5(654  4-  3652)*)  4-  (9653  +  5655)2)  4-  5(65*  4-  6)1)2)  4-  (96)*  4- 15653)2)  4-  (56)3  4-  8654)2)  4-  5(65*  4-  365i)2)] 
=  5[5(654  4-  8653)2)  4-  5(6)3  +  55*)*)  +  (965*  4-  565i )*)  4-  5(6)*  4-  3653)2)  4-  (96)3  4-  5654)**  4-  (965*  4-  1565i  )2)) 
x»l»i  =  (9654  +  15652)*;  +  9(6)3  +  553)*)  +  (56)*  4-  96)1  )*)  4-  (96,*,  4- 1561*3)*)  +  (SS^,  +  96)4)2)  4-  5(6)*  4-  36), )2) 

*^7«  —  (86)4  4- 156)3)*)  +  (56)3  4-  965s)2j  4-  9(6)*  4-  6)1  )*)  4-  5(6)*  4-  8653)2)  4-  (96)^  4-  56)4)2)  4-  (96)*  4-  156)i  )2) 

’*i««  ~  5(654  +  3633)*!  4-  (5615  4"  96)3)23  4-  (963*  4-  5641  )23  4-  (96**  4-  15653)*i  +  8(653  +  5),)*)  4-  (96)*  4-  156)i  )27 

u^9i  =  (36)4  4"  5632)*!  4"  3(6*3  4"  6)5)23  4-  3(6**  4-  6*1  )23  4-  (36**  4-  56)3)2)  4-  8(6)3  +  5*, )*)  4"  (36)*  4-  56)i )27 

t«£)i  =  (86)4  4- 156)3)*;  +  (86)7  4- 1563*3)2)  +  (56)*  4-  964*1)2)  4- 15(64*3  4-  6)*)*)  4-  6(6)3  +  363*1)2)  4-  6(6)7  +  36)4)2) 

Wilt  =  5(674  +  36*3)*!  4-  5(6)7  4-  36*3)23  4"  (96*8  4-  56*1  )24  4- 15(6*3  +  5)*)*)  +  (86*3  4-  156)i  )*)  4-  (86)7  4- 15634)2) 

••te*  ®  (5674  +  56*3)21  4"  (36i7  4-  5633)23  +  3(6)*  4"  6,1  )*,  4-  5(6)3  4"  6**)*)  ■+  (36*3  +  56*1)2*  4-  (86*7  +  56*4)2) 

**^31  —  5(637  4-  655)21  4-  (367*  4-  56)3)23  4"  5(634  +  5)1  )*3  4-  (36*4  +  66)3 )*)  4-  (86*7  4-  56)5)2*  4"  5(6)*  4-  6*1  )2) 

**44«  =  5(6)4  +  573)*;  +  5(6)*  4-  6)3)2)  4-  (36)7  4-  56)4)2)  4-  (36)*  4-  66)1)2)  4-  6(6)7  4-  6)5)2)  +  (36)  4-  66)1)2) 

ti^si  =  15(634  +  553)21  4-  (9673  4"  156)3)23  4"  5(6*7  +  8631)2)  4-  (56*3  +  86)4)24  4-  (96)7  4"  156)5)2)  +  6(6*5  +  86)1)27 
xx>2»i  —  15(634  4"  653)*!  4- 5(6*3  4- 36)3)23  4- (96*7  4- 156*1  )*)  4- (96*3  4- 56)4)*)  4- 5(6)7  4"  36)5)25  4- (96*5  4-  156)i)*7 
*^7«  —  5(634  +  653)21  4"  (3673  4"  56)3)*)  4"  (36*7  +  56)1)2)  4-  3(653  4"  6)4)*)  4"  (86)7  4-  66)5)2)  4"  (36*  4*  56)1)2} 


14 


and  the  geometric  constants  a*,,  ,  and  6|l^„  are  given  by 


C‘1 

= 

6m 

/JS. 

=; 

Cm 

-^5 

«m 

— 

— 

''fn 

Cm 

_ 

“m 

= 

6m 

ai? 

= 

Cm 

6m 

Cm 

6m 

Cm 


SU 

sL 

6i„ 


The  values  for  zj  for  i  =  1,2,3,. ..,8  are  given  by 


6nCm  ~  6n»Cn 
Cfi^m  —  Cm^ 


i 

^2 

-4- 

24 

26 

2? 

1 

-1 

1 

-1 

-1 

1 

-1 

2 

1 

-1 

-1 

1 

-1 

1 

-1 

3 

1 

1 

-1 

-1 

1 

-1 

-1 

4 

-1 

-1 

1 

1 

1 

-1 

-1 

5 

-1 

1 

-1 

1 

-1 

-1 

1 

6 

1 

-1 

1 

-1 

-1 

-1 

1 

7 

1 

1 

1 

1 

1 

1 

1 

8 

-1 

-1 

-1 

-1 

1 

1 

1 

(90) 


(91) 


The  stiffness  coefficients  due  to  compatible  displacements  are  then  given  explicitly  as 


K^..=gtiguj  ;  j  =  1,2,3,. ...24;  i=  1,2,3 . t;  ib  =  l,2,3 . 6N,  (92) 

In  computing  the  stiffness  contributions  of  the  incompatible  displacements,  the  selected  nonconforming 
displacement  modes  are  identical  to  those  presented  in  [12]  and  are  given  by. 


where 


■  Ml 

M2 

Ms 

0 

0  0 

0 

0 

0 

Ai 

A2 

0 

0 

0 

Ml 

A/2  A/3 

0 

0 

0 

i 

0 

0 

0 

0 

0  0 

Ml 

M2 

M3 

A9 

Ml  =  e-M-f2r}-f3(: 
M2  =  fi^-M-Av-ZsC 
Ms  =  C*-M-/8»7-/9C 


The  derivation  of  the  incompatible  modes  and  all  constrmts  are  contained  in  Appendix  III. 


(93) 


(94) 


The  components  of  equation  (31),  r,-;- ,  where  t  =  1, 2, 3, ...,  Ne  and  j  =  1, 2, 3  are  given  by 


where 


rij  =  diiGlj 

rij+3  =  di20? 

>*ij+6  =  dizQfj 

fi+Ar.j  =  daiQjj 
’‘i+N.j+3  = 

’^i+N.J+e  =  dj30,.y 


ri+2N..j 

^i+2N,j+3 

^i+2N.j+6 

n+3N,J 

ri+3N,J+3 

ri+3N.j+6 


=  d310(, 
=  ^3205 
=  d330j 
=  0 
=  d440^ 
=  d440.?; 


^i+4N„} 

''»+4W.j+3 

’■•+4W.J+8 

'•t+SN.J 

’■•+5W.J+3 

ri+3N,j+6 


=  d550?j 

=  0 
=  d550l 
=  deM 
=  ^6601 
=  0 


e‘ 

>=1 


(95) 


(96) 
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«*i  =  8[/,  (6*,  +  3^1^  +  /2(i  *6  +  3^*3)  +  MSti  +  )  +  2(6^^  +  i*3)]/3 

<I>U  =  8(A(«5‘4  +  3^1,)  +  hiSU  +  3if3)  +  M6L  +  3«,*,)  +  2(i|i  +  tfS4)]/3 

^?3  =  8[/7(65^  +  3^^)  +  ^(6,^  +  3^‘3)  +  ^(6Js  +  3^5l)  +  2(^5^  +  ^*6)]^ 

For  i  =  2, 3, 4 

=  [16«;,^i  +  40(u;,^4/i  +  tyfj/j  +  ti>f6/3)]/45 

=  [16u;/'j  +  40(ti>,\/4  +  u/fs/s  +  ti», ^6/6)1/45 

+  AQ(wfj7  +  tufj/s  +  ti;f<j/9)i/45 

For  I  =  5, 6, 7, 8, 9, 10 

^?l  =  -[V-?1+/i^?4  +  /2V'?5-/3V'?6]/135 

=  -[0?2  +  M?4  +  /5V'?6-/6<<'?6l/135 

«?3  =  -[V'f3  +  /4V'f4  +  /50?5-/6^^f6]/135 

The  constants  Vij  and  are  given  by 

V>Ji  =  48o\ti>f7  +  16(a'2ti;fio  -  4»"ni  +  “>*12  “  “5“'a4  “  “6“'*i6  “  “V*"*!?) 

V»Jj  =  48aju;?8  +  16(aVu»?j4  -  -  a‘5«;?27  +  o;,ti>,^9  -  +  aitwfjs) 

^.•3  =  48a^u;{'9  +  16(0*5 tuf28  -  «6«'*30  +  «7«'‘32  +  “iW'iM  -  «2«"*95  -  “iu'fss) 

=  24o‘>?i+  8(a*7u;?i6-a>f,o  +  o*5u;f„)  + 40(0*8  tu, ^18 +  a^u>fi5-Qriu;f,3) 

V»?5  =  24oit«f2-  8(o*7u;?27-a>f,3  +  o^u;f2e)-40(a^wfj2  +  o*5tu?25-a*i‘",-2o) 

^>*6  =  24o^ti;,^3+  8(oiiu?3<-oJu;f35-o*7«>f38)  +  40(o\ii;f29-a2“'*3i-“4“'.‘33) 


(97) 


(98) 


(99) 


(100) 


*"fi  =  ^(9^f»  +  156*3)  +  503(637  +  6*4)  +  5/14(673  +^45)  +  /1s(3647  +  5695)  +  ^4(3675  +  5633) 

*®?2  —  5/Jj(6*i  +  6*4)  +  ^3(96*4  +  156*1 )  +  504(6*7  +  6*j)  +  ^5(36*4  +  56*j)  +  07(5637  +  56*3) 

w*3  =  5/33(6*4  +  6f7)  +  501(673  +  6**)  +  ^4(96*8  +  156*2)  +  03(56*7  +  56*4)  +  07(56*3  +  56*2) 

w*4  =  3/93(634  +  6*2)  +  /Js(^T4  +  5632)  +  /^U^ST  +  3633)  +  /9i(6?3  +6^4)+  03(6*7  +  6^) 

®?9  =  ^(^4T  +  36fj)  +  503(6*3  +  6^5)  +  ^4(6*8  +  36J3)  +  ^s(6*5  +  6^7)  +  07(6*3  +  6m) 

=  ^(^78  +  36*1 )  +  03(^«7  +  36*4)  +  501(637  +6*i)  +  03(6*3  +  6*1 )  +  07(6*7  +  6m) 

—  503(6*7  +  6m)  +  501(6*2  +  6m  )  +  dU56*7  +  56m)  +  01(56*4  +  56m)  +  ®^2(^34  +  6*3) 
w*3  =  507(6*7  +  633)  +  501(6m  +  6?i )  +  1^1(3687  +  56m)  +  ^03(6*4  +  ^41 )  +  ^(56ji  +  36*4 ) 
u7*9  =  507(6*7  +  644)  +  503(6^3  +  6*1)  +  90'4(633  +  6ei )  +  01(5637  +  5634)  +  ^(36*5  +  56*1 ) 
wfio  =  /8»(36*r  +  56m)  502(56*4  +  5642)  +  507(6*7  +  6*3)  +  501(56^  +  56*3)  +  15/9^6*3  +  6*4) 

w?xi  =  ^4(36*7  +  5634)  +  502(56*3  +  5633)  +  507(6*7  +  6^4)  +  503(56*3  +  56^)  +  15^4(6*5  +  6*3) 

tp*i3  =  3^4(637  +  6m)  +  ^7(36*7  +  56m)  +  ^03(6*3  +  ^44)  +  503(56*4  +  56J3)  +  901(6*2  +  6m) 

*®?J3  =  fi»(6*7  +  ^m)  +  07(6*7  +  36m)  +  502(6*3  +  6*4)  +  501(6*4  +  56*2)  +  501(6^4  +  6m) 

*®?14  =*  3/85(6*7  +  6m)  +  ^7(36*7  +  56m)  +  995(6*3  +  644)  +  994(6*4  +  9*3)  +  395(36^4  +  56^)  (101) 

w?is  —  01(637  +  9*4)  +  95(9*7  +  39m)  +  395(9*3  +  94s)  +  501(6*3  +  6*3)  +  394(9*4  +  39*3) 

wfi4  =  ^(56*3  +  56*3)  +  94(39*4  +  56^3)  +  95(39*7  +  59m)  +  501(633  +  9*3)  +  595(937  +  9*4)  +  595(934  +  9^3) 

^i\7  =  3/84(633  +  9*3)  +  03(56*3  +  5633)  +  95(39*7  +  56*4)  +  395(9*3  +  633)  -f  395(9*7  +  9*4)  +  95(36*4  +  56*3) 

W*18  *  03(633  +  6^3)  +  03(6*3  +  5633)  +  95(9*7  +  39^4)  +  01(633  +  673)  +  03(637  +  9*4)  +  95(9*4  +  39m) 

ti>*i9  =  395(9*7  +  944)  + /3»(3647  +  56*4)  +  993(9*1  +  9m)  +  395(3644  +  56*1)  +  995(9*4  +94i) 


wfso  =  95(9*7  +  633)  +  01(637  +  56*3)  +  395(644  +  6*1 )  +  395(644  +  56*1 )  +  395(9*4  +  941) 

t»j3i  =  394(917  +  6*4)  +  03(56*7  +  5614)  +  995(646  +  9*3)  +  995(9i6  +  9*3)  +  395(39*4  +  5613) 

®>22  =  95(9i7  +  9*4)  +  94(647  +  39^4)  +  395(6*3  +  644)  +  395(6*4  +  9*3)  +  395(9m  +  36*3) 

'•’ias  “  03(5637  +  59*4)  +  393(39*4  +  56*1 )  +  595(9*7  +  633)  +  395(59*i  +  3644)  +  15/85(9*4  +  637 ) 

Wj34  *  394(641  +  6*4)  +  95(3647  +  56*4)  +  95(56*1  +  36m)  +  395(9*7  +  644)  +  395(6*1  +6*4)  + 0i (56^4  ■+56*7) 
®?2s  =  95(9*4  +  9*3)  +  03(6*3  +  39*3)  +  95(9*7  +  39j‘4)  +  95(9*3  +  9m)  +  95(9*7  +  9*4)  +  95(9m  +  39f3) 
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and  the  geometric  constants  a*,,  0^,  and  S^„  are  given  in  (90). 

The  computation  of  the  inner  product  given  in  (32)  results  in  a  9  x  9  Z  matrix  which  must  be  inverted 
to  give  the  coefficients  aij .  Although  Z  is  symmetric  and  specific  explicit  inversion  schemes  can  be  applied 
to  minimize  computations,  this  matrix  order  is  perhaps  at  the  limit  in  which  closed-form  expressions  for  the 
inverse  may  be  succintly  expressed  and  a  numerical  scheme  may  be  preferred.  The  stiffness  contributions 
due  to  the  incompatible  contributions,  A'a,  are  given  by 

ft  A,,  —  ~9ni^m^tm^km9kj 


where 


»=  1,2,3 . 6Ne  \  j  =  1,2,3,...,!  ;  n,fc=l,2,3 . 24;  s,m  =  1,2,3, ...9 


The  complete  element  stiffness  matrix  is  therefore  given  by  the  sum 

A'  =  A'^  +  A'a 


Tke  E8NL  Element 

The  E8NL  element  is  formulated  using  stress  expansions  defined  in  natural  coordinates.  Incompatible  modes 
are  introduced  to  complete  the  quadratic  bending  terms  and  enforce  orthogonality  between  the  stress  and 
incompatible  strain  fields  variationally.  Two  versions  are  presented  to  demonstrate  the  effect  of  different 
stress  expansions  on  element  behavior.  ESNLi  is  based  on  incomplete  quadratic  expansions  for  the  inplane 
stress  components  with  complete  linear  expansions  for  the  shear  stresses  which  is  identical  to  the  FEl  ele¬ 
ment  presented  in  Reference  [10].  The  selected  incompatible  displacement  modes  used  in  E8NL  are  identical 
to  those  used  in  the  E8PL  element. 

The  assumed  stress  field  is  given  by 


r  = 


[ri] 


[ri 


(rij 


(raj 


[rz] 


[rz] 


(102) 


where 


(103) 


The  natural  stresses  are  transformed  to  physical  coordinates  using  a  contravariant  transformation  based  on 
centroidal  Jacobians  and  premultiplied  by  the  distributing  matrix  to  yield  the  transformed  stress  held  as 

P  =  D-‘r,r  (104) 

Without  attempting  a  simplification  of  the  modes  in  (104)  through  linear  combinations  of  the  unknown 
expansion  coefficients,  the  orthonormalizing  process  yields  stress  modes  of  the  general  form 

Pil  +Phi  +P*3»1  +Ph^n  +Pi6'fC  +Pf7<€ 

P*a  +Pi9i  +P.*io»?  +  P.*uC  +  P.’i2^»H-P,*i3»K  +  P*m« 
p»  _  ,  PilS  +  PilS^  +  Pa?*?  +  PilsC  +  Pil9^^  +  PiJO’jC  +  Pi’ziC^ 

*  Pin  +  Pi23^  +  PUaV  PiSsC 
P?36  +  PiZT^  +  Pii»^  +  PinC 
.  Piao  +  Pr31^  +  Phi’}  +  PiwC 

where  the  coefficients  are  obtained  numerically.  The  components  =  1,2, 3,...,  33  and  n  = 

1,2,3,  ...,8  result  in  the  following  lineau  algebraic  expressions 

7 

lf(m)(3n-2)  =  X^(‘^nP(m,*)  +  <^2lP(m,t+7)  +  ‘^3lP(m.*+14))®n*  + 

t=l 
4 

^(d55P(m.t+25)®nt  +  ‘^66P(fn,*+29)®nt) 
k=l 
7 

ff(m)(3n~l)  =  +  ‘^22P(m.li+7)  + ‘^32P(m,t+14))®nt  +  (106) 

i;=l 
4 

^(d44P(m,t+21)®nl:  +  ‘^66P(m,fc+29)®n*) 

fcssl 

7 

Sl(m)(3a)  =  X^(<^13P(m,fc)  +  <^23P(*m,t+7)  +  ‘/33P(m.t+l4))©n*  + 

t=l 
4 

^(d44P(m,t+2l)^nt  +  ‘^55P(m,t+25)®nt) 
t=l 


where 

e*„,  =  -(30j.,+^U  +  K9  +  <^^io)/3  eu  =  -(3^U  +  <>j.2o)/27 

e^2  =  -(3^1.2  + <13 +  ^’ni5)/9  ei,e  =  -(3^1.6 +  <i«)/27 

®n3  —  ~(3^n3  +  <ll  +<I6)/0  ©1,7  —  "(3^1,7  +  ^nl9)/27 

Qn4  =  ~(3^n4  +  ^nl2  +  ^nI4)/0 

and 


(107) 


<1  ~  ^"^32  —  ^5.^!3  ~  ^12^7 

<2  =  *"(^34  “  ^52)  +  ^2^13  ~  ^3^*12  +  ^5  ~  ^7^*14 

<3  —  *"^62  +  *2^32  +  ^5(^43  +  ^16.)  ~  ^6^12  ~  *7^42 
<4  —  *r^36  +  *3  <2  +  ^5^53  +  ^6  ^13  ~  ^7  (^16  +  ^m) 
<5  =  *?(«j4  +  <2)  +  +  <3)  +  +  ^’17)  - 

*3  ^42  ~  ^4<2  ~  ^6  ^1.4 

^nS  ~  *6  (^43  “  ^52)  ~  77(^46  +  ^72)  +  3^5  (^73  +  ^ss)  + 
*2  +  *3  <2  +  7"  ^32 

^n7  =  3:" (^37  +  tfjg)  +  23(^34  —  ijg)  -  2” (5*17  +  ^J^)  + 

*2  +  76^15 

<8  =  *"^54  +  72^15  —  73  ^14 
<9  ~  ^2^2  ~  75  ^46  ~  ^6^42 


<10  —  736^  +  Zg  {33  —  zy^jg 
<11  =  ^1^74  +  ^2^17  ~  74^14 
<12  =  ^"^57  ~  ^3^17  +  z"6}5 
<13  =  32^72  —  Zg  j^2  +  ^3^47 
^ril4  =  3:g5g2  +  Zgifg  —  Zg{72 
<15  =  ^3<7  ~  ^7^57  +  74^53 
^nlS  =  ^4  ^36  .+  ^6  ^73  ~  ^7^6 
<17  =  2(2” ^5g  +  Z3  jgg  +  Zg  igg) 

<18  =  ^2^57  +  ^3^74  +  34^45+ 
<19  =  3:2^76  +  Z"6g4  +  23^47 
^n20  —  ~  ^3  ^76  +  *4  ^56  —  Zg  ^57 


(108) 


18 


wumpmm 


The  values  for  z"  are  given  in  (91). 

The  components  rmn.  with  m  =  1,2,3 . 33  and  n  =  1,2,3,  are  given  by 

•"m.n  "b  *^2l6mn2  "b  ^3l6rnn3  "b  ^S5^mn2  d"  ^66^mn3 

»‘m.n+3  =  <fl20mnl  +  <^220mn2  +  <^326mn3  +  +  <^66^mn3 

•"min+S  ~  *fl3  0mnl  +  <^236mn2  +  ‘^33©mn3  +  +  <^55*mn2 


(109) 


where  the  fdlowing  constants  are  defined 


"ml* 

Qm2*  = 
ei„3*  = 


i>Uik 

^m2k 

V'mS* 

V>m4/b 

V>mSt 

^m6k 

V'mT* 

^m9k 

'l^mllk 


-^inlk  +  M\nAk  +  M\nkk  +  Mm6k  ^ml* 
-  ^Uik  +  MU*k  +  A  '^mS*  +  A  V-jnet  *m2k 
^m3*  +  MmAk  +  M\nkk  +  Mm6k  ^Lsk 


~^m7k  +  A^mlOt  +  AV’mll*  +  AV’ml2t 

-^mSt  +  A^mlOt  +  AV-^IU  +  f«^ml2k  (110) 
^m&k  +  h^m  10*  +  AV’Inll*  +  AV’mi2t 


3p(m.7*-6)A4 

3p(m,7*-8)-^5 

3P(m,7*-6)-^6 

3p(m,7*-5)A4 

3P(m.7t-8)-^2 

3p(m.7*-8)-^‘l 

3P(m.4*+18)'^i 

3p(m,8*+18)-^5 

3p(m,4*+18)<^6 

3p(m,8*+l9)'^4 

3P(m,4*+20)‘^5 
3p(  m,  8t+18)M 


+  P(m.7t-2)‘^22  +P(m 
+  P(m.7t-2)-^23  +P(m 
+  P(m,74-4)-^16  +P(m 
+  P(m.7*-3)-^*17  +P(m 
+  P(m.7*-5)-^8  +P(m 
+  P(m,7t-5)'^7  +  P(m 
+  P(m,4*+20)'^U+P( 

+  P(m.4*+20)A*i7+  P( 

+  P{m,4*+19)A‘i2+  P( 

+  P(m,4t+21)'^‘l7+  P( 

+  P(m,4*+19)-^8  +P( 

+  P{m,4t+19)‘^7  +P( 


7t-4)-^14  +P(m,7t)A26  +  P(m,7t-3)-^‘l9  +  P(mJk-S)K 

7fc-4)'^‘lS  +P(Fn, 74-1)^24  +  P(m,7t-3)'^20  +  P(m.7t-5)-^10 
7*)-^27  +  P(m,7*-l)-^25  +  P(m,7t-3)'^21  +  PCm^t-S)-^!! 

7t-4)'^*13  +  2p(m,7t-2)-^14  +  2P(in.7t)A\9  +3p(m,7t_6)A3 

,7t-3)'^18  +  2p(m,7t-l)A20+2p(m,7t-2)A‘io  +  3p(m,7*-4)A5 
,7*-4)-^‘i2  +  2p(m,7*)A\i  +3p(m_7t_3)A6  +2p(mj*_  i)A',5 


.4t+21)-^19+  P(m.4*+19)‘^9 
,4t+21)-^23+  P(fn,4*+I9)<^ll 
,4*+21)'^24+  P(m,4t+20)'^18 
,4t+20)'^13'l'  3P(m,4*+18)-^3 
,4t4-21)'^‘l8+  3P(m,4t+18)‘'2 

.4*4-20)  A|2+  3p(m,4*4-2I)A6 


(111) 


A’l 

= 

8(«56+3«|2)/9 

A' 

0 

=  16(47  +  45) 

A'i9 

16(47  +  46)/9 

A* 

= 

8(4  +  3«‘i3)/9 

A' 

1 

=  16(47  +  44) 

A2O 

= 

16(48 +  43)/9 

= 

8(«i4  +  3f^,)/9 

A' 

2 

=  8(46-1-342) 

A?i 

(14446  +  22b6y/Ab 

Ai 

= 

8(^34  +  ^52)/9 

A' 

3 

=  8(^74  -1-  342) 

A22 

= 

(4844  +  8042)/45 

A’s 

= 

8(«;3+«i8)/9 

A' 

4 

=  16(42  +  44) 

A23 

= 

(4847  +  8045)/45 

(112) 

Ai 

= 

8(«i6  +  4)/9 

A 

5 

=  (14446 -»-22543)/5 

A24 

= 

(4846  4-8043)/45 

A? 

= 

8(«b  +  3^i4)/9 

A 

6 

=  16(48  +  42) 

Am 

(48«^6  +  806i2)/45 

As 

= 

8(«i7  +  3«i5)/9 

A 

7 

=  8(47-l-3i$6) 

A26 

= 

(4847-I-80«^)/45 

Ai 

(14444  +  22542)/45 

A' 

8 

=  8(46  +  343) 

A27 

= 

(4847  +  80«54)/45 

and 

^mn  ~  —  f>m 

Cjt 

t 

^mn  ~  CfiOni  —  i 

6L 

k  “ 

—  nm4 

(113) 

The  componente  Oij  are  then  obtained  through  the  inversion  of  the  symmetric  9  x  9  Z  matrix.  The  el¬ 
ement  stiffness  matrix  is  given  explicitly  by 


where 


^ij  ~  "I"  —  9ni9nj  ~  9ni7'n$^tm^km9k}  i  ^ji  — 

1=1, 2, 3 . 24;  i  =  1,2,3 i;  n,fc  =  1,2,3 33;  s,m  =  1,2,3, ..., 9 
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A  variation  of  the  above  element  formulation,  designated  ESNLo,  is  presented  to  show  the  simplification 
possible  if  balanced  expansions  are  used  for  the  stress  components.  The  assumed  stress  field  is  given  by 


r  = 


[ri] 

[r,] 

[r.l 

[ri] 


(114) 


where 

Ti  =  [1.^.  (115) 

The  natural  stresses  are  transformed  to  physical  coordinates  using  a  contravariant  transformation  based  on 
centroidal  Jacobians  and  premultiplied  by  the  distributing  matrix  to  yield  the  transformed  stress  field  as 

P  =  D-‘r„r  =  TT  (116) 

Through  a  linear  combination  as  defined  by 

0  =  diag[T,  f .  f .  f ,  (117) 

the  P  matrix  can  be  reduced  to  the  uncoupled  form  of  stress  modes  given  in  (114).  Performing  the  orthonor- 
matixation  of  the  fundamental  modes  yields 


P  —  [pI>P2>P3<P4iP5>P6>P7l 


(118) 


where,  for  i  =  1,2,3,  ...,7,  the  general  form  of  each  mode  is  given  by 

Pi  =  Pii  +  +  Pis';  +  PUC  +  Ph^n  +  Pis  nC  +  PirCi  ( 1 1 9) 


and 

Pifc  =  0  k  >  i 


Explicit  expressions  for  the  coefficients  pj^  are  presented  in  Appendix  III. 

The  components  gmn  w’th  m  =  1,2, 3,...,  7  and  n  =  1,2, 3,...,  8  result  in  the  following  linear  algebraic 
expressions 


7  7  7 


^(m)(3n— 2) 

9(m)l3n-l)  =  ^  ^nPmk^nk 

9(m)(3n)  =  Y^dl3Pmk^lk 

tsl 

t=i 

7 

7 

7 

9(m-|-7)(3n-2) 

fl'(m+7M3n-I)  =  ^•^liPmk^nk 

P(m+7)(3n)  =  ^<i23Pmk^lk 

k=l 

k=l 

k=l 

7 

7 

7 

9(m-t-14)(3n-3) 

i;(m+H)(3n-l)  =  Y^d23Pmk^nk 

ff(m+M)(3«)  =  y^<f33PmtQnt(^20) 

k=l 

•9 

t=l 

S(m+2lK3n-2) 

=  0.0 

4 

ff(m+21)(3n-l)  =  ^‘^44Pm*©n* 

7 

ff(m+21)(3n)  =  ^**Pmk^lk 

7 

k=l 

ksl 

7 

tf(m+28)(3n-2) 

=  y^<^5SPm.t0L 

l!  =  l 

'/(m+28)(3n-l)  =  00 

P{m+28)(3n)  =  ^  <^SSPm,k^hk 
fcsl 

i;(m4-3S)(3n-2) 

k=l 

7 

fl(m+aS)(3n-l)  =  ^66Pm.k^hk 

k=l 

i;(m+35)(3n)  =  0.0 

20 


where  0^^,  and  rj*  are  defined  above  in  equations  (107),  (108),  and  (91),  respectively. 


The  components  r^n 

with 

m  =  1,2,3,. 

..,7  and  n  = 

1,2,3,  are  given  by 

’’m.n 

duQLn 

**m,n+3 

di2^mn  ’’m.B+S 

di3eL 

^m+7,n 

dl2&mn 

^m+7,n+3 

d22^mn  ^m+7,n+6 

= 

d23e^n 

dlsOLn 

^ni+M,fi+3 

= 

d23^mn  ^m+14,fi+6 

= 

d33ei„ 

^m+21,n 

Z= 

0.0 

l^m+21  ,n^3 

= 

^44^mn  ^m+21,n+6 

= 

d44ei,„ 

^m+38,n 

= 

d55ei„ 

^fn+28,fi+3 

2= 

0.0  »‘m+38,n+6 

desefl,„ 

^m+35,n 

dee^mn 

^fn+35,n+3 

= 

deeOjnn  *‘m+3S,n+6 

= 

0.0 

where  the  following  constants  are  defined 

©ml  = 

+  /2^m5  +  /al/’mS 

©m2  = 

+  Mins  +  Mine 

0i„3  = 

+  MmA 

+  MinS  +  Mine 

V'ml  =  +  Pm5  A22  +  Pm3-^H  +  Pm7^26  +  Pm4‘^19  +  Pm2'^? 

V'mJ  =  ^pJljjA's  +  PmsAja  +  PmsAjs  +  PmS'^M  +  Pm4-^20  +  Pm2‘'*10 

^m3  —  3pmi  Afi  +  PmsA'is  +  Pm7^27  +  Pm6'^25  +  Pm4'^21  +  Pm2-^*11 

^•^4  =  3Pm2A^  +  Pm^A'jy  +  PmsA’ia  +  2pJ),5A‘,4+  2pJ),7A*i9  +  3pJ),i  A3 

V’ms  =  3pJ),,  A2  +  PmjAg  Pm4X\a  +  2pJ),6A2o+  2pJ),5A\o  +  3p^3A^ 

^in6  =  3p;;.iA‘  +  p;;.2A‘7  +  p;;.3A\2  +  2p;;,7A* ,  ^  3p;;„A^  +2p;;.6A‘6 

and  the  constants  At  and  are  defined  above  in  equations  (112)  and  (113),  respectively. 
The  element  stiffness  matrix,  A',  is  given  by 


(121) 


(122) 


(123) 


where 


Axj  —  1^4, j  +  A’\,j  —  <Jni9nj  ~  9ni^n$<^tm^km9kj  i  Aji  —  I\ij 

i  =:  1,2,3,...,24  ;  j  =  1,2,3, ...,i  ;  n,t  =  1,2,3, ..., 42  ;  s,m  =  1,2,3,  ...,9 


8  Numerical  Demonstrations 

In  order  to  validate  the  performance  of  the  explicitly  derived  4  and  8-node  elements,  several  standard  bench¬ 
mark  problems  are  analyzed.  Solutions  to  plane  stress  problems  are  presented  in  Figures  5  and  6  while 
Figure  7  depicts  solutions  to  a  solid  cantilevered  beam  problem. 
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Figure  5.  5-eiement  cantilevered  beam  subjected  to  (1)  pure  bending  and  (2)  end  shear. 
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Figure  6.  Clamped  circular  beam  under  end  shear  loading. 
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Figure  7.  Solid  cantilevered  beam  under  end  moment  loading. 


9  Conclusion 

The  aim  of  introducing  special  stress  field  transformations  has  been  to  maximize  the  efficiency  of  hybrid 
element  formulations  by  allowing  explicit  forms  of  element  stiffness  matrices  to  be  derived.  By  fully  ex¬ 
ploiting  the  freedom  in  selecting  and  transforming  independent  stress  fields  in  the  hybrid  stress  technique, 
the  computational  cost  associated  with  numerical  matrix  integration  can  be  eliminated  and  inversions  can 
be  reduced  substantially.  In  the  extention  to  higher-order  element  formulations  such  as  the  8-node  plane 
and  20-node  solid  elements,  without  introducing  incompatible  displacement  modes,  matrix  inversions  are 
eliminated  entirely.  The  approach  has  been  demonstrated  by  deriving  explicit  linear  algebraic  forms  for 
the  stiffness  matrices  of  selected  4-node  quadrilateral  and  8-node  hexahedral  elements.  The  computational 
advantage  over  purely  displacement-based  element  formulations  is  clearly  evident  and  the  method  outlined 
in  this  study  should  be  expected  to  find  genered  application  in  various  hybrid/mixed  methods  to  minimize 
computations  in  determining  element  stiffness  coefficents. 
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APPENDIX  I 


DETERMINATION  OF  DISTRIBUTING  MATRICES 

In  computing  the  2-D  distributing  matrix,  D,  the  inversion  of  the  material  compliance  matrix,  5,  is  not 
required  as  the  D  matrix  may  be  directly  related  to  the  material  stiffness  matrix  C  as 

where  the  C  and  D  matrix  are  given  for  an  orthotropic  material  as 


C  = 

The  eigenvalues  are  computed  as 
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—  (  \/^22  2ci  1C22  +  ^*^12  "I"  ^1 1  "f"  ^22  "t"  1  )/2 
1^2  ~  (  \/^22  “  2ci  1 C22  +  4c^2  "I"  •'1 1  +  ^22  +  Cii)/2 
fa  =  C33 


yielding  the  matrix  as 

The  Q  matrix  is  defined  as 
where  the  eigenvectors  are  given  by 


Q  = 


f  C.2  ] 

(p2  —  C22 

f  ®  ) 

<  Fi  -  Cl, 

If 
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II 
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1  0  J 

i  0  J 

,  C33  J 

and  normalized  as 


=  Ni^i 


The  computation  of  the  3-D  distributing  matrix  using  the  symmetric  C  matrix  for  an  orthotropic  material 
IS  performed  as  follows. 
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The  eigenvalues  are  obtained  as 
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where 


<1  =  (r+vV  +  rS)‘/3 

h  =  (r  - 

q  =  b/Z-a-/9  ;  r  =  (ab  -  3c)/6  -  a^/27 

a  =  -(C33  +  C22  +  Cii) 

l>  =  {(C32  +  Cil)C33  -  C23  “  <^13  +  C11C22  —  cfj) 

C  =  ((C11C22  +  Ci2)C33  +  C11C53  -  2Ci2Ci3C23  +  cf3C22) 

yielding  the  matrix  as 


The  Q  matrix  is  given  by 


Q  -  [*ll<t2|*3|*4|*5|*6l 


where  the  eigenvectors  are  determined  by 
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and  are  normalized  to  yield 

Ni = 

=  Ni^i 


In  the  case  repeated  eigenvalues  corresponding  to  a  coupled  submatrix  partition,  the  eigenvectors  asso¬ 
ciated  with  the  degenerate  eigenvalues  are  discarded  and  are  replaced  by  vectors  orthonormalized  to  the 
independent  eigenvectors  using  the  Gram-Schmidt  procedure. 

APPENDIX  11 

DETERMINATION  OF  STRESS  MODE  COMPONENTS 
USING  THE  GRAM-SCHMIDT  PROCEDURE  ■ 

An  elaboration  of  (22)  yields  an  automated  procedure  for  the  Gram-Schmidt  process  and  is  presented  to 
indicate  the  numerical  simplicity  of  developing  orthonormal  stress  polynomials.  A  generic  algorithm  may 
be  developed  by  first  defining  variables  representing  the  initial  stresses,  intermediate  combinations  and  the 
final  orthonormalized  stresses  given  by  pijk,  p'ijk  and  respectively,  where  t  indicates  the  mode  number, 
j  denotes  the  stress  component  and  k  represents  the  coefficient  in  the  polynomial  expansion  for  the 
component.  Next,  the  basic  scalar  integrals  involving  the  Jacobian  determinant  and  the  polynomial  powers 
arising  from  the  inner  product  defined  in  (21)  are  computed  and  assigned  to  a  variable,  accounting  for 
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the  order  of  expaosion  and  the  arreuigement  of  terms.  For  example,  in  the  2-D  case,  if  the  expansion  for 
each  stress  component  is  given  by 

Pil  +Pi2*  +  Pi3tf 

Pi  =  ■ 


PiA  +  PiS*  +  Pi6» 
Pt7  +  PiS*  +  PiiV 


1 


where  z  and  y  are  expressed  in  physical  or  natural  coordinates,  the  scalar  integrals  required  are  given  by 

,1  ,1  f  1  I  r  Ai  Aj  As 

^  =  /  /  (UK  *  >[l. *. jiUi'fy  =  Aj  A4  As 

I  y  J  L  ^3  . 


The  computation  of  scalar  integrals  is  detailed  in  Appendix  IV. 

With  ^  defined,  the  orthonormaiizing  procedure  is  given  by  the  following  sequence  of  operations  which  may 
be  performed  symbolically  or  numerically 


i-i 

0  Ptmn  —  P»mfi  ~  ^_^(Prt«P»l:j  )Prmn 

r=l 


**)  P$mn  ~  PsmniPtkiPtkj^ij) 


As  shown,  the  nature  of  the  Gram-Schmidt  procedure  is  such  that  in  developing  an  orthonormal  set,  each 
mode  is  computed  in  an  ’evolutionary’  fashion  through  the  process  of  orthogonalizing  the  mode  to  the 
preceding  t  —  1  modes.  When  performed  symbolically,  this  leads  to  increasingly  cumbersome  expressions 
for  coefficients  as  a  function  of  the  number  of  modes  and  degree  of  stress  coupling  and  is  most  efficiently 
implemented  as  a  numerical  procedure.  Explicit  expressions  are  presented  below  as  examples  for  select  2-D 
elements  to  indicate  the  unwieldly  nature  of  symbolic  representation. 

For  the  E4PQ  element,  the  coefficients  in  the  orthondrmalized  stress  modes  are  given  by 
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Normalizing  constants  are  given  by 
Ni  = 

=  [a?(A5-Al/Ai)  +  ai(A,-A2/Ai)]-‘/’ 
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and  material  constants  are  given  by 
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For  the  E4PR  and  E4PL  elements,  the  coefficients,  Pij,  of  the  stress  modes  are  given  by 
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For  the  EI4NL  element,  the  stress  mode  coefficients,  p,j,  we  given  by 
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For  the  E8NL2  element,  the  coefficients,  pij ,  of  the  stress  modes  are  given  by 
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where 
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Ns  =  [RlXi  +  R^\3  +  R^X9-^Pi\io  +  Xia  +  2(RsR6\2  +  RsR7X3-RsPsX4  +  (Rs  +  RiR7)X3- 
R6f^X^  +  R6X14  —  RtPsXs  +  /Z7Ai2  —  /’sAll)]"’^" 

Ne  —  [•^8^1  +  R10X9  +  iZjjAio  +  P9  Ajs  +  Aig  +  2{R3R9X2  +  R3R10X3  4-  R6R11X4— 

(RsPi  +  R9Rio)Xs  +  (Ra  +  ^lo^iOAe  +  /Zg/ZuAr  —  R9PaXi4  +  {R9  —  iiiiP9)Aii  — 

R10P9X13  +  ilioAiT  +  AiiAis  -  P9A23)]”‘^* 

JV7  =  [PijAi  +  i2i3A8  +  P14A9  4-  RisXiQ  +  Pi4Xi9  +  PfgAis  +  A20  +  2(/ii2iil3A2  +  Pl2Pl4A34- 
P12P15A4  +  (PiaPie  +  RizR’i*)^^  +  (RuRii  —  Pi2Ri4)Xs  +  (^12  +  Ri3Ris)X7+ 
AiaAieAu  +  (RisRie  +  Pm  -  Pi3Pm)Ah  +  P13A16  +  PmPi6Ai2  -  PmPmAi7- 
PisPmAis  +  P15A13  —  Pi6PmA22  +  P16A23  —  PmA2i]~*^* 


APPENDIX  III 

INCOMPATIBLE  DISPLACEMENT  MODES 

The  constants  computed  for  the  element  incompatible  displetcement  modes  are  based  on  the  approach  pre¬ 
sented  in  [12]  for  identically  satisfying  the  convergence  criterion 

LMdv  =  0 

for  any  arbitrary  element  configuration. 


Tht  El  PL  and  E4NL  elements 


For  the  E4PL  element,  the  various  constants  used  in  the  definition  of  the  incompatible  displacements  are 
given  by 


0163  —  ai6| 
atbj  —  02^1 


f2  = 


0363  —  0362 
0162  ~  ®2^1 
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3ai63  +  0263(72  —  fx) 
36163  +  6263(72  -  fx) 

3(63  +  61 72  -  6270 
0i037i  ~  30102  O2O372 

6i037i  —  36ia2  —  62O372 


3(027i  —  Oi72  —  03) 

016371  —  30261  —  026372 
616371  —  36261  —  626372 
3oi02  +  O2O372  ~  Oi037i 
361O2  +  62O372  -  6i037i 


For  the  E4NL  element,  7i  wid  72  ore  the  same  as  given  above.  The  constants  hj  are  slightly  different  and 
are  given  by 

hi  =  8(63  +  6i72  —  627i)/3  he  =—80372/9 

h2  =  8(363  —  6371  )/9  /»7  =  86372/9 

63  =  863/2/9  he  =  — 8(36i  +  6372)/9 

he  =  8(02/1  —  03  —  oi72)/9  hg  =—803/1/9 

hs  =  8(03/1  —  3o2)/9  hio  =  8(301+0372)79 


The  E8PL  and  E8NL  elements 

A  comment  on  the  approach  in  [12]  of  forming  incompatible  modes  has  been  made  in  [17]  which  criti¬ 
cizes  the  algebraic  complexity  in  the  3-D  case.  However,  with  careful  manipulation,  the  procedure  is  quite 
tractable  and  is  presented  in  full  detail  below. 

The  basic  incompatible  displacements  are  selected  as 

«A  =  K^^^C-] 

and  ’virtual  parameters’  are  taken  as 

“A  =  K,«?,C] 

The  elements  pij  of  P,  are  obtained  from  the  integration  of 


ml  <11  <12  <13 

<21  <22  <23  d^dridC 

.  <31  <32  <33 


where  <y  =  Adj(J)ij  and  which  yields 

Pii  =  8(045  +  3a23)/3:  P21 
P12  =  8(064  +  3o3i)/3;  P22 

Pi3  =  8(056  +  3oi2)/3;  ,P23 

where 


The  inverse  of  P,  is  given  by 


8(^45  +  3^23)/3;  P31 
8(/?64  +  3/?3i)/3;  P32 
H0i6  +  3/?12)/3;  P33 


8(645  +  3623)73 
8(654  +  3631)73 
8(656  +  36i2)/3 


oy 

=  6,Cj  -  bjCi 

0ii 

—  CjOj  —  Cjdx 

h 

—  Oibj  —  Ojbi 

0-1  -  J_ 

*  ~\P.\ 


P22P33 

-  P23P32 

PI3P32 

-  P12P33 

PI2P23 

-  P13P22 

Pll 

P12 

Pl3 

P23P31 

-  P21P33 

PuP33 

-  P13P31 

P13P2I 

-  PIIP23 

= 

P2I 

P22 

P23 

P21P32 

-  P22P3I 

P12P3I 

—  PI1P32 

P11P22 

-P12P21 

. 

P32 

P33 

where 


jft  I  =  Pll(P22P33  -  P23P32)  -  Pl2(P2lP33  -  P23P31 )  +  Pl3(P2lP32  -  P22P3l) 

The  elements  py  of  Px  are  obtained  from  the  integration  of 

/I  /•!  /■!  r  ^<11  ¥12  C<13  Pll  P12  Pl3 

/  /  ^<21  ¥22  <<23  d^dqdC  =  P21  ft2  P23 

xJ-lJ-X  y^t3i  T,<32  <<33 


Pll  P12  Pl3 

iff 

P21  P22  P23 

*  I  f 

P31  P32  P33 
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which  yields 

Pu  =  16(«25  +  a«)/3  Ppi  =  16(^5  +  043)/^  =  16(^35  +  ^431/3 

P}7  ~  +  0'34)/3  P22  =  16(^1  +  034)/^  PSJ  =  16(561  +  634)/^ 

Pi3  —  16(ai6  +  flf5j)/3  P23  =  16(^16  +  /?5i)/3  P33  =  16(6i6  +  582)/3 

The  computation  of  the  final  form  of  is  given  by 


such  that 
where 

and 


UA  = 

[Mi.A/o.Ms] 

/ « \  1 

■  Ml 

M2 

Mz 

0 

- 

0 

0 

0 

0 

6 

V 

II 

e 

II 

0 

0 

0 

Ml 

M  2 

Mz 

0 

0 

0 

V- A 

0 

0 

0 

0 

0 

0 

Ml 

M2 

M3 

Ml  =  e-M-hf}-f3c 
M2  =  r}^-  /4^  -  hn  -  feC 
Mz  =  e-M-hn-M 


The  constants  fi  in  the  definition  of  the  incompatible  modes  are  given  by 


/i  =  PiaPai +PI2P21 +  PuPii;  h  =  P23P31  +  P22P21  +  P21P11;  /a  =  ^^Pai  +  W2P21  +  WiPn 

/4  =  PT3P32  +  P12P22  +  PUPl2i  /*  =  P23P32  +  P22P22  +  WlPl2>  /s  =  rt3P32  +  1^2P22  +  P3lPl2 

/t  =  Pi3P33  +  PnP23  +  PHP13;  /a  =  P23P33  +  P22P23  +  P2lPt3;  /»  =  PMP33  +  P32P23  +  P31Pi3 


APPENDIX  IV 

COMPUTATION  OF  BASIC  SCALAR  INTEGRALS 


In  the  Gram-Schmidt  procedure,  the  inner  product  defined  by  (21)  requires  the  evaluation  of  scalar  in¬ 
tegrals  in  generating  the  orthonormal  basis  vectors.  For  the  2-D  elements,  using  complete  linear  stress 
expansions  in  physical  coordinates,  the  required  integrals  consist  of  polynomial  orders  up  to  quadratic  order 
defined  as  ^  ^ 

[Ai,  A2,A3,  A<,A5,A6]  =  IX  \J\[i,x,y,x^,y^,xy]didrj 

In  computing  the  above  integrals,  the  determinant  of  the  Jacobian  is  expressed  as 

|J|  =  Jo  +  Ji^  +  J211 


where 

Jo  —  0162  “<*2^1  i  Ji  —  0163  —  0361  ;  J2  —  0362  —  0263 
which  under  mapping  to  isoparametric  coordinates  by 


X  =  oi4  -I-  0217  -I-  03^q 

y  =  6i^  -P  6217  +  63^*7 


31 


yield  closed  form  expressions  given  by 
Aj  =  AJq 

Aj  =  4(Jiai  +  J2a2)/3 

A3  =  4(Ji6i  +7262)73 

A4  =  [4Jo(3oi  +  3^2  +  **3)  +  3(7in2<*3  +  720103)1/9 

As  =  [47o(36i  +  36j  +  63)  +  8(716263  +  726163)1/9 

As  =  4(63(7003  +  7i03  +  720i)  +  62(80270  +  037i)  +  6i(3oi7o  +  0372)1/9 


For  complete  linear  stresses  expressed  in  natural  coordinates  the  basic  integrals  are 

[Ai,A2,A3,As,A5,A6l  =  y  j  17|[l,^,fj,f»j,^*,»7*ldC7i/ 

which  yield  simply 

Ai  =  47o  Ail  =  0 

A2  =  47i/3  As  =  47o/3 

A3  =  472/3  As  =  47o/3 

In  the  3-D  case,  the  determinant  of  the  Jacobian  is  given  by 


|7|  =  ri  +  r2^ +  r3Jj  +  r4C  +  r5j;4  +  rsC^  +  rrCJJ  +  rs^* +  r9»;* +  rio<* +  riiC»j{ +  ri2»;f^+ 

+  '’isC'/*  +  +  rntiC  +  +  '■2o^C* 

where 


ri 

= 

»*6 

= 

V’li  +  'Pai  + 

*•11 

=  ^pU 

ne 

=  '^as 

ra 

= 

Vai  +  ¥>12 

= 

'^23  4"  <^*34  +  (^62 

’’12 

=  p\a 

n? 

=  Pai 

ra 

s 

v’m  +  ‘Pij 

r» 

= 

‘PU 

»’13 

-  Psi 

ns 

=  ^54 

U 

= 

<pIi  +  via 

*•9 

= 

pIz 

**14 

=  ^42 

ns 

=  P46 

ra 

= 

V’42  +  ^14  +  ¥»12 

»'10 

Pas 

ns 

=  V?26 

no 

=  Pes 

and 

=  Oj(oy6fc  -  akbj)  +  bi{bjCt  -  b^Cj)  +  Cj(c;Ot  -  ctaj) 

With  assumed  stresses  of  quadratic  order,  the  basic  integrals  required  in  computing  the  orthonormalized 
stress  modes  are  defined  up  to  quartic  order  by 


[Ai,A2,A3,  As,  As,A6,  Ar,  As,  Aol  =  /:/:/:  |7|[1,  z,  y,  z,  zz,  yz,  zy,  z*,  y'ld^drjdC 

[Aio.Air,  Ai2,Ai3,Ai4,Ais,  Aisl  =  y  j  J  \J\[z^,xyz,xy^,xz^,yx^,yz^,zx^]d^dr)d(; 


[Ai7,-'i8,''i9,A20,A2i,A22,A23l  =  y  J  J  |7([zy^, z^ y^, z^z, y^ ZZ®, yz^lrffdi/dC 

[A24,  A2S,  A28,  A27,  A28.  A29,  A30I  =  LLL  |71[yr®,  zz®,  zy®,  z®y® ,  y®z®,  z®z® ,  zyz®ldf  dqdC 
[A311 A32,  A33,A34,A35l  =  j  J  J  \J\[xy^z,x^yz,x*,y*,z*]d^dT)dC 


where  the  isoparametric  mapping  is  given  by 


=  oi^  +  02q  +  03C  +  a4^q  +  os^C  +  osqC  +  07^»?C 
=  6i^  +  62T/  +  63C  +  64^*7  +  65^C  +  66qC  +  67^»7C 

=  Cl^  +  CnT]  +  C3C  +  C4^T7  +  Cs^C  +  CsT/C  + 


X 

y 

z 
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Closed  form  expressions  for  the  integrals  may  be  succintly  given  for  bwer  orders,  however,  for  arbitrary  orders, 
the  integrals  of  the  scalar  functions  are  most  expediently  evaluated  using  a  standard  Gaussian  integration 
scheme. 

N.  N, 

A„  =  5;  ^  2  Aijt  |i(«. .  Vi .  <:k)\x'Ui,  Vi .  Ct  )y'(6 .  Vi ,  Ck)z'{ii ,  Vi .  Ck) 

tsl  j  =  l  tsl 

For  stresses  assumed  in  natural  coordinates,  the  basic  integrals  required  in  ccwnputing  orthogonal  spanning 
modes  are  of  simple  form  and  can  be  evaluated  explicitly  for  any  arbitrary  powers  in  terms  of  coefficients, 
r,-,  of  the  Jacobian  determinant.  For  the  E8NL  element,  23  basic  integrals  are  defined  by 


•  1  ' 

Ai 

Ag 

A3 

A4 

As 

As 

A7 

A2 

As 

As 

A7 

Ai4 

An 

Ais 

f*  f* 

V 

A3 

As 

Ag 

Ag 

Ai2 

Ai7 

All 

[|J| 

C 

’  [i.^.v.C^v.vCXd^^^vdC  = 

A4 

A7 

Ag 

Aio 

All 

Ais 

Ai3 

y_i  y_i 

^v 

As 

Ai4 

Ai2 

All 

A 18 

A22 

A23 

vC 

As 

All 

Ai7 

Ais 

A22 

Ai9 

A21 

«  J 

A7 

Ais 

All 

Ai3 

A23 

A21 

Ago 

and  evaluate  to 


Ai  =  8(3ri  +  rg  +  rg  +  rio)/3 

Ag  =  8(3r2  +  ri4  +  ri6)/9 

A3  =  8(3r3  +  ri2  +  rn)/9 

A4  =  8(3r4  +  ri3  +  ri5)/9 

As  =  8(3r5  +  r2o)/27 

As  =  8(3r7  +  ri8)/27 

A7  =  8(3rs  +  ri9)/27 

As  —  8(15ri  +  9r8  +  Srg  +  5rio)/45 

Ag  =  8(15ri  +  brg  +  9rg  +  5rio)/45 

\io  =  8(15ri  +  Srs  +  5rg +  9rio)/45 

\ii  =  8rii/27 

\i2  —  8(15r2  +  9ri4  +  5ri6)/135 


Ai3  ■=  8(15r2  +  5ri4  +  9ri6)/135 

Ai4  =  8(15r3  +  9ri2  +  5ri7)/135 

Ai5  =  8(15r3  +  5ri2  +  9ri7)/135 

Ais  =  8(15r4  +  9rj3  +  5ri5)/135 

At7  =  8(15r4  +  5ri3  +  9ri5)/135 

Ai8  —  8(15ri  +  9r8  +  9r9  +  5rio)/135 
Ajg  =  8(15ri  +  Srg  +  9r9  +  9rio)/135 
A20  =  8(15ri  +  9r8  +  Sfg  +  9rio)/135 
A21  =  8(5r5  +  3r2o)/135 
A22  —  8(5r6  +  3rig)/135 
A23  —  8(5r7  +  3ri8)/135 
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